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1. 


INTRODUCTION 


The  U.S.  Navy  needs  accurate  gravity  field  data  for 
navigation  and  missile  guidance  applications.  The  scientific 
community  also  uses  marine  gravity  data  in  studies  of  litho¬ 
spheric  structure  and  global  tectonics.  Most  of  the  necessary 
gravity  field  information  cannot  be  measured  directly  and  must 
be  derived  computationally  from  large  volumes  of  existing  meas¬ 
urement  data.  Such  data  may  come  from  diverse  sources  --  for 
example,  satellite  altimeter  measurements  as  well  as  conven¬ 
tional  ship  survey  gravity  data. 

Existing  data  bases  of  measured  gravity  field  data 
are  expected  to  expand  considerably  in  the  next  few  years  as, 
for  example,  data  from  the  gravity  gradiometer  survey  system 
become  available.  Large  quantities  of  data  continue  to  be 
received  from  the  GEOSAT  radar  altimeter  mission.  Emphasis 
must  be  placed  now  on  the  development  and  evaluation  of  compu¬ 
tational  techniques  to  estimate  geodetic  quantities  of  interest 
to  the  Navy  and  the  scientific  community.  These  techniques 
will  have  to  exploit,  in  an  accurate  and  computationally  effi¬ 
cient  manner,  an  ever  expanding  data  base  including  many  types 
of  measurements,  each  with  its  own  noise  and  spatial  distribu¬ 
tion  characteristics.  The  necessary  computational  techniques 
must  also  provide  realistic  accuracy  estimates  for  the  gravity 
quantities  being  determined. 

The  International  Association  of  Geodesy  (IAG)  has 
recently  established  a  Special  Study  Group  to  investigate  algo¬ 
rithms  for  local  gravity  field  determination.  This  group  has 
be  in  examining  the  accuracy  and  computational  efficiency  of  a 
number  of  geodetic  algorithms,  using  a  standard  set  of  gravity 
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field  data.  Algorithms  are  being  examined  and  evaluated  for 
efficiency,  stability,  accuracy,  and  applicability,  as  well  as 
their  limitations.  Study  group  participants  have  conducted 
numerical  tests  of  the  computation  of  gravity  anomalies,  deflec¬ 
tions  of  the  vertical,  geoid  undulations,,  and  gravity  gradients. 
The  group  has  also  considered  issues  of  error  estimation  and 
incorporation  of  topographic,  geophysical,  and  geological  infor¬ 
mation  into  the  estimation.  TASC  has  participated  in  the  work 
of  this  group,  primarily  in  the  area  of  applying  the  GEOFAST 
algorithm  as  a  tool  for  the  efficient  determination  of  gravity 
field  quantities. 

The  concept  of  geodetic  algorithm,  illustrated  in 
Fig.  1-1,  is  central  to  the  work  described  in  this  report. 

This  term  includes  methods  for  estimating  and  computing  gravity 
field  quantities  of  interest  on  the  basis  of  available  measured 
data.  The  gravity  field  quantities  of  interest  (geoid  undula¬ 
tion,  gravity  anomaly,  deflection  of  the  vertical,  gr^ity 
gradients,  etc)  may  not  be  directly  measurable  or  not  measur¬ 
able  conveniently  and  economically.  This  definition  of  geo¬ 
detic  algorithm  is  presented  in  more  detail  in  Section  2.1  of 
the  report;  an  introduction  to  the  major  categories  of  geodetic 
algorithms  is  given  in  Section  2.2.  Appendix  A  is  a  detailed 
review  of  these  algorithms. 

The  objective  of  this  report  is  to  present  the  results 
of  a  study  to  determine  the  most  appropriate  techniques  and 
algorithms  for  estimating  geodetic  quantities  from  a  variety 
of  observation  types.  The  investigations  reported  here  have 
included  both  theoretical  and  computational  studies.  Figure  1-2 
summarizes  these  investigations. 

This  report  includes  seven  chapters.  Chapter  2  is 
devoted  to  an  examination  of  geodetic  algorithms.  In  addition 
to  an  overview  of  the  concept  of  a  geodetic  algorithm,  the 
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chapter  introduces  the  classical  integral  evaluation  methods, 
as  well  as  least-squares  collocation,  GEOFAST,  and  Fourier 
transform  methods.  Details  of  all  of  these  techniques  are 
reviewed  in  Appendix  A. 
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Chapter  3  presents  a  summary  of  the  methods 
IAG  Study  Group  investigations.  Comparative 
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of  various  algorithms  and  techniques  are  presented  in  text  and 
tabular  form.  This  chapter  summarizes  the  key  results  on  the 
relative  accuracy  and  efficiency  of  the  methods  considered  by 
the  IAG  Study  Group. 

Chapter  4  summarizes  various  aspects  of  the  treatment 
of  topographic  or  terrain  effects,  especially  important  in 
regions  of  high  relief  such  as  the  White  Sands  study  area.  In 
addition  to  reviewing  the  relevant  efforts  and  results  of  the 
IAG  Study  Group,  Chap^r  4  includes  an  overview  of  original 
investigations  cor-4  .  ._ed  by  TASC ,  with  a  more  detailed  dis¬ 
cussion  in  Append  C. 

Chapter  5  describes  the  essentials  of  a  new  derivation, 
developed  by  TASC,  of  a  Fast  Fourier  Transform  (FFT)  algorithm 
for  implementing  terrain  correction.  The  detailed  derivation 
of  this  approach  is  documented  in  Appendix  B. 

In  Chapter  6,  results  are  presented  for  a  series  of 
numerical  experiments  involving  applications  of  the  GEOFAST 
algorithm.  The  two  major  phases  of  this  study  involved  compu¬ 
tations  using: 

•  Synthetic  marine  gravity  data  derived 
from  the  high-degree-and-order  part  of  a 
global  geopotential  model 

•  A  set  of  synthetic  high-frequency  local 
gravity  data  constructed  on  the  basis  of 
the  statistical  characteristics  of  local 
gravity  truth  data. 

The  principal  conclusion  of  this  study  is  that  modern  geodetic 
algorithms  can,  with  appropriate  data  input,  resolve  deflec¬ 
tions  of  the  vertical  to  0.2  arc  sec  for  the  marine  data  and 
to  0.5  arc  sec  for  the  local  gravity  data. 
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Chapter  7  presents  TASC's  conclusions  and  recommenda¬ 
tions  resulting  from  the  study  eff-irt. 

In  addition  to  the  summaries ,  discussions,  and  evalu¬ 
ations  of  the  work  of  the  IAG  Study  Group,  and  the  presenta¬ 
tion  of  appropriate  supporting  material  from  other  relevant 
published  scientific  literature,  this  report  documents  several 
original  contributions  by  TASC  to  the  geodetic  algorithm  study 
effort.  These  contributions  include: 

•  Analysis  of  GEOFAST  algorithm  performance 
(Chapter  6) 

•  New  and  improved  derivations  of  FFT  ter¬ 
rain  correction  algorithms  (Chapter  5 
and  Appendix  B) 

•  Studies  of  the  prediction  of  local  deflec¬ 
tions  of  the  vertical  from  terrain  height 
data  alone  (Chapter  4  and  Appendix  C). 
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OVERVIEW  OF  GEODETIC  ALGORITHMS 


This  chapter  introduces  the  concept  of  geodetic  algo¬ 
rithm  and  includes  technical  discussions  of  several  main  cate¬ 
gories  of  such  algorithms.  Section  2.1  defines  geodetic 
algorithms  and  suggests  a  variety  of  classifications  based  on 
criteria  which  include  type  and  distribution  of  available  meas¬ 
urement  data,  quantities  to  be  computed,  and  the  statistical 
or  deterministic  orientation  of  the  computational  approach. 
Section  2.2  introduces  a  number  of  specific  algorithms  which 
are  presented  in  more  detail  in  Appendix  A.  These  include 
classical  integral  evaluation  methods,  least-squares  colloca¬ 
tion,  TASC's  GEOFAST  algorithm,  and  other  Fourier  transform 
and  frequency  domain  approaches. 


2.1  GEODETIC  ALGORITHMS 

The  term  geodetic  algorithm  is  used  extensively 
throughout  this  report  to  refer  to  methods  of  computing  geo¬ 
detic  quantities  of  interest  from  available  measured  data. 
These  algorithms  can  be  classified  and  categorized  in  a  number 
of  ways : 

•  The  specific  geodetic  quantitj^  to  be 
computed  or  estimated 

•  The  general  approach  (deterministic  or 
statistical ) 

•  The  type  or  types  of  measurement  data  to 
be  used 
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•  The  geographical  distribution  of  the 
points  at  which  the  quantities  are  to  be 
computed  --  for  example,  a  small  number 
of  irregularly  spaced  points  as  opposed 
to  a  large  number  of  points  forming  a 
regular  grid 

•  The  distribution  of  the  points  at  which 
data  measurements  are  assumed  to  be 
available . 

In  addition,  algorithms  can  be  distinguished  on  the  basis  of 
specific  techniques  used  for  the  numerical  evaluation  (exact 
or  approximate)  of  the  underlying  theoretical  formulations. 

A  brief  discussion  of  the  various  bases  of  classifi¬ 
cation  follows.  It  is  intended  to  be  an  introduction  to  the 
large  variety  of  specific  algorithms  that  are  treated  in  the 
rest  of  this  report. 

Geodetic  quantity  to  be  computed  -  Historically,  the 
earliest  examples  of  geodetic  algorithms  were  developed  to 
compute  the  undulation  of  the  geoid  or  the  components  of  the 
deflection  of  the  vertical  at  a  specified  location.  In  the 
first  case,  the  undulation  is  a  geodetic  quantity  that  cannot 
be  measured  directly,  even  if  the  specified  location  is  easily 
accessible.  In  the  second  case,  while  the  deflection  of  the 
vertical  can  be  measured  directly  at  an  accessible  point  using 
astrogeodetic  techniques,  there  are  substantial  economic  and 
logistical  considerations  that  may  make  indirect  determination 
more  attractive.  Also,  many  locations  are  inaccessible  in  a 
practical  sense  to  the  equipment  required  for  direct  astro¬ 
nomic  measurement  of  the  deflection. 

Other  geodetic  data  types  which  have  become  important 
through  more  recent  applications  include: 
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•  Components  of  the  gravity  disturbance 
vector  at  altitude 

•  Components  of  the  gravity  gradient  tensor 

•  Mean  values  of  the  gravity  anomaly,  typi¬ 
cally  for  regions  within  which  direct 
gravity  measurements  are  scarce  or  un¬ 
available  . 


There  is  a  general  unifying  concept  that  facilitates 
the  consideration  of  the  numerous  geodetic  quantities  that  may 
be  the  subject  of  computation  or  estimation  algorithms.  This 
point  of  view  regards  the  disturbing  potential  as  the  fundamen¬ 
tal  geodetic  quantity  to  be  determined  from  whatever  combination 
of  observational  data  may  be  available;  other  geodetic  quanti¬ 
ties  are  computed,  in  principle,  by  applying  linear  operators 
of  various  kinds  to  the  disturbing  potential  (Fig.  2.1-1). 
Thus,  for  example,  Bruns's  formula: 

N  =  T/y  (2.1-1) 


where 

N  is  the  undulation  of  the  geoid 
Y  is  the  normal  gravity 
T  is  the  disturbing  potential 

is  the  simplest  example  of  this  concept.  Similarly,  the  gravity 
disturbance  vector  is  expressed  as: 

5  =  grad  T  (2.1-2) 

the  scalar  gravity  disturbance  as: 

6g  =  -  |I  (2.1-3) 


2-3 


THE  ANALYTIC  SCIENCES  CORPORATION 


Figure  2.1-1  Relationships  Among  Geodetic  Quantities 

where  the  partial  derivative  is  taken  along  the  normal  to  the 
geoid  (plumb  line);  and  the  gravity  anomaly  as  (exact  form): 


Ag 


+  I  ?1  T 

3n  y  3n 


(2.1-4) 


or,  in  the  usual  spherical  approximation: 


3T  2 
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where  the  partial  derivative  is  taken  in  the  radial  direction, 
and  R  is  a  mean  radius  of  the  earth.  This  concept  is  treated 
in  more  detail  in  Appendix  A. 2. 3,  where  such  operators,  illus¬ 
trated  in  Fig.  2.1-1,  are  summarized  in  Table  A. 2-1. 

Deterministic  and  statistical  approaches  -  Classical 
geodetic  algorithms  assume  that  the  measured  quantity  is  known 
without  error  everywhere  on  the  globe.  For  example,  the  Stokes 
integral  for  determining  undulation  or  the  Vening  Meinesz  inte¬ 
gral  for  deflections  of  the  vertical  both  require,  in  theory, 
error- free  knowledge  of  the  gravity  anomaly  everywhere .  Actual 
evaluations  of  these  algorithms  involve  finite  approximations 
to  infinite  processes  like  integration,  as  well  as  estimation 
(interpolation  or  extrapolation)  of  unmeasured  gravity  anomaly 
data  in  terms  of  available  measurements. 

Modern  statistical  approaches  to  the  formulation  of 
geodetic  algorithms  attempt  to  deal  directly  with  issues  in¬ 
volving: 

•  The  discrete  and  finite  nature  of  the 
measured  data 

•  The  statistical  nature  of  measurement 
error  in  the  data 

•  The  formulation  of  best  methods  based  on 
global  statistical  properties  of  measured 
data 

•  The  development  of  methods  using  a  com¬ 
bination  of  data  types  in  an  optimal 
manner. 

Collocation  methods,  discussed  in  Appendix  A. 2,  are  an  example 
of  statistically  oriented  geodetic  algorithms. 
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Measurement  data  types  -  Classical  geodetic  algorithms 
are  based  on  the  use  of  point  gravity  anomaly  data,  since  this 
is  the  type  of  geodetic  measurement  that  can  most  readily  be 
made  directly  and  reasonably  economically.  Modern  algorithm 
development  is  marked  by  the  availability  of  a  wide  variety  of 
other  data  types,  such  as: 

•  Geopotential  field  coefficients  recovered 
from  satellite  orbit  studies 

•  Satellite  altimeter  estimates  of  undula¬ 
tion  over  the  oceans 

•  Astronomic  measurements  of  deflection  of 
the  vertical 

•  Mean  gravity  anomaly  data  of  specified 
block  size 

•  Direct  measurements  of  components  of  the 
gravity  gradient  tensor 

•  Gravity  field  quantities  recovered  from 
other  high-accuracy  position-determination 
technologies,  such  as  inertial  and  GPS. 

Algorithm  development  is  also  influenced  by  the  desirability 
of  using  mixed  data  types  as  input  for  the  determination  of 
other  geodetic  quantities  of  interest. 

Distribution  of  computation  points  -  Classical  algo¬ 
rithms  were  designed  to  produce  results  for  a  single  point  at 
a  time,  with  the  bulk  of  the  computations  having  to  be  redone 
entirely  to  obtain  a  result  for  another  point,  even  if  it  is 
nearby.  (However,  some  exceptions,  such  as  the  method  of  astro- 
gravimetric  interpolation,  can  be  noted.)  A  characteristic  of 
many  modem  algorithms  is  the  ability  to  produce  results  es¬ 
sentially  simultaneously  for  a  number  of  points,  which  may  be 
distributed  irregularly  or  located  on  a  uniform  rectangular 
grid. 
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Distribution  of  data  points  -  Classical  algorithms 
assume,  at  least  in  principle,  the  availability  of  data  every¬ 
where  on  the  surface  of  the  earth.  Modern  algorithms  relax 
this  unrealistic  assumption  in  various  ways.  Collocation  meth¬ 
ods  (see  Appendix  A. 2),  for  example,  are  based  on  the  optimal 
use  of  a  finite  number  of  irregularly  distributed  observations 
of  various  kinds,  corresponding  closely  and  without  idealization 
to  the  realities  of  data  measurements  in  the  field.  The  GEOFAST 
algorithm  (see  Appendix  A. 3),  on  the  other  hand,  operates  on  a 
regular  rectangular  grid  of  data  points  and  outputs  values  on 
a  regular  grid. 


2.2  INTRODUCTION  TO  SPECIFIC  ALGORITHMS 

The  algorithms  documented  in  detail  in  Appendix  A 
include  classical  integral  evaluation  methods,  least-squares 
collocation,  TASC's  GEOFAST  algorithm,  and  other  Fourier  trans¬ 
form  and  frequency  domain  approaches.  A  brief  overview  and 
discussion  of  each  of  these  methods  is  given  in  the  present 
section . 


Integral  evaluation  methods  -  The  classical  integral 
evaluation  methods  of  Stokes  and  Vening  Meinesz  determine  the 
undulation  of  the  geoid  and  the  deflection  of  the  vertical  at 
a  specific  point  from  global  free-air  gravity  anomaly  data. 
More  recent  modifications  of  these  methods,  due  to  Molodenskii 
and  Pellinen,  make  explicit'  use  of  terrain  height  data  and 
refer  to  the  physical  surface  of  the  earth  rather  than  to  the 
geoid.  Both  the  classical  methods  and  more  recent  developments 
are  discussed  in  Appendix  A.l. 

The  classical  methods  may  eventually  be  supplanted  by 
modern  approaches  based  on  the  use  of  gravity  disturbances 
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rather  than  gravity  anomalies,  associated  with  precise  local 
positioning  techniques  such  as  the  use  of  the  Global  Position¬ 
ing  System  (GPS).  In  addition,  the  computational  approach  of 
numerical  evaluation  of  integrals  taken  over  the  earth's  sur¬ 
face  may  largely  give  way  in  practice  to  the  collocation  and 
frequency-domain  techniques  outlined  below.  Nonetheless,  the 
classical  integral  evaluation  methods  are  important  for  at 
least  three  reasons: 


•  They  are  the  theoretical  bcsis  for  more 
recent  developments 

•  They  are  still  widely  used  in  practice 
in  a  number  of  areas  (e.g.,  Ref.  6;  see 
also  Ref.  49) 

•  They  were  considered  to  be  competitive 
in  the  White  Sands  Study  (Ref.  21). 


Least-squares  collocation  -  This  modern  technique  for 
estimating  geodetic  quantities  from  observed  data  is  discussed 
in  Appendix  A. 2.  Among  the  advantages  it  offers,  as  compared 
with  classical  integral  evaluation  methods,  are  the  following: 


•  The  use  of  a  variety  of  measured  data 
types  as  input,  and  the  simultaneous 
estimation  of  several  gravity  field 
quantities 

•  The  use  of  irregularly  spaced  observa¬ 
tion  points 

•  The  incorporation  of  statistical  informa¬ 
tion  about  errors  in  the  input  data, 
with  corresponding  results  for  the  qual¬ 
ity  of  the  estimated  geodetic  quantities 

•  The  estimation  or  improvement  of  model 
parameters  along  with  the  determination 

of  the  geodetic  quantities  to  be  estimated. 
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An  important  characteristic  of  the  collocation  approach  is 
that  least-squares  estimation  is  applied  to  data  from  which 
systematic  and  long-wavelength  components  have  been  removed. 
These  residual  data  and  the  quantities  estimated  from  them 
correspond,  therefore ,  . to  a  more  nearly  stationary  and  iso¬ 
tropic  residual  gravity  field.  As  a  final  step,  the  appro¬ 
priate  systematic  and  long-wavelength  components  are  restored. 

The  GEOFAST  algorithm  -  TASC's  GEOFAST  algorithm 

(Ref.  16)  was  developed  as  a  computationally  efficient  method 

for  solving  the  minimum-variance  estimation  equations  that 

arise  from  the  application  of  least- squares  collocation  to  the 

determination  of  geodetic  quantities.  It  replaces  the  direct 

matrix  solution  of  the  estimation  equations,  a  process  whose 

3 

computational  cost  is  of  the  order  of  N  (where  N  represents 
the  number  of  data  points),  by  a  sequence  of  operations  for 
which  the  overall  computational  effort  is  of  order  N  log2  N. 
The  GEOFAST  algorithm  is  described  in  Appendix  A. 3;  Chapter  6 
of  this  report  presents  the  results  of  two  sets  of  numerical 
experiments  --  one  using  marine  gravity  data  in  an  ocean  trench 
area,  the  other  using  land  gravity  data  in  the  White  Sands 
area . 


Other  Fourier  Transform  Methods  -  Methods  based  on  a 
two-dimensional  Fourier  transform  lead  to  efficient  computa¬ 
tional  algorithms  that  estimate  geodetic  quantities  in  a  vari¬ 
ety  of  ways.  Examples  include  the  evaluation  of  the  Stokes 
integral  for  undulation  of  the  geoid,  the  Vening  Meinesz  inte¬ 
gral  for  deflection  of  the  vertical,  the  Poisson  integral  for 
upward  continuation,  and  the  convolution  integrals  defining 
terrain  corrections.  Appendix  A. 4  presents  the  theory  of 
Fourier  transform  methods  and  describes  specific  applications 
to  convolutions  and  to  the  evaluation  of  Stokes  and  Vening 
Meinesz  integrals. 
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3.  SUMMARY  AND  REVIEV  OF  IAG 

STUDY  GROUP  METHODS 


This  chapter  concerns  two  reports  on  the  application 
of  geodetic  algorithms  in  mountainous  terrain.  The  first 
(Ref.  20;  also  to  be  referred  to  as  the  "1983  Report")  was 
produced  by  members  of  the  IAG  Special  Study  Group  4.70  on 
"Gravity  Field  Approximation  Techniques."  Several  sets  of 
investigators  contributed  descriptions  of  the  application  of  a 
variety  of  algorithms.  Two  main  tasks  were  undertaken:  the 
prediction  of  gravity  anomalies  and  the  prediction  of  deflec¬ 
tions  of  the  vertical  from  gravimetry  and  gridded  topographic 
heights  in  the  White  Sands,  NM  area.  The  second  report 
(Ref.  21;  "1985  Report")  was  written  jointly  by  a  group  of 
five  investigators  gathered  at  the  University  of  Calgary  and 
concentrates  solely  on  the  prediction  of  deflections  of  the 
vertical.  Four  techniques  are  described  and  compared  in  detail 

The  chapter  is  organized  into  three  sections.  Sec¬ 
tion  3.1  describes  briefly  the  White  Sands  data  sets;  illustra¬ 
tions  and  more  complete  descriptions  are  included  in  Chapter  4 
and  Appendix  C,  which  describe  TASC ' s  further  analysis  of  these 
data.  The  gravity  anomaly  prediction  algorithms  from  the  1983 
Report  are  summarized  in  Section  3.2.  The  prediction  of  deflec 
tions  of  the  vertical  Js  described  in  Section  3.3;  it  incorpo¬ 
rates  material  from  both  reports. 


3.1  WHITE  SANDS  DATA 

The  IAG  study  group  used  four  main  sets  of  data  from 
the  White  Sands  area.  A  brief  description  of  each,  along  with 
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its  geographic  location  and  number  of  data  values,  is  given  in 
Table  3.1-1.  The  gravity  anomaly  data  are  divided  into  two 
groups  to  facilitate  the  evaluation  of  algorithms  for  predict¬ 
ing  gravity  anomalies;  the  second  set  contains  truth  data. 

The  locations  of  the  points  in  data  sets  1  and  2  that  fall 
within  the  3  deg  by  3  deg  region  between  32  and  35  deg  North 
and  105  and  108  deg  West  are  shown  in  Fig.  C.2-1.  Their  values 
were  gridded  to  form  the  color  image  in  Fig.  C.2-2.  The  grid- 
ded  heights  (data  set  3)  are  illustrated  as  a  color  image  in 
Fig.  C.2-3,  for  the  region  between  32  and  35  deg  North  and  105 
and  108  deg  West.  The  vector  plot  of  Fig.  C.5-1  shows  the 
locations  and  values  of  the  astrogeodetic  deflections  of 
data  set  4.  Note  that  in  the  plots  of  these  data  shown  in 
Fig.  4  (p.  7)  of  the  1983  Report  and  Fig.  1.6  (p.  11)  of  the 
1985  Report  the  sign  of  the  east  deflection  n  is  reversed  from 
the  usual  convention  (e.g.,  Ref.  1). 


TABLE  3.1-1 
WHITE  SANDS  DATA  SETS 


TYPE 

LOCATION 

NUMBER  OF  DATA 

Point  free-air 

31-36.4  deg  N 

11,277 

gravity  anomaly 

104-109.4  deg  W 

(with  height  and  location) 

Point  free-air 

33-34  deg  N 

102 

gravity  anomaly 

106-107  deg  W 

(with  height  and  location) 

Point  heights  on  a 

31-35  deg  N 

172,320 

regular  grid,  30  arc  sec 

105-108  deg  W 

spacing 

Pairs  of  astrogeodetic 

mainly 

deflections  of  the 

32-34  deg  N 

441 

vertical 

106-107  deg  W 
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3.2  GRAVITY  ANOMALY  PREDICTION 

The  study  group's  strategy  to  evaluate  algorithms  for 
predicting  gravity  anomalies  was  to  use  data  set  1  (free-air 
gravity)  supplemented  by  data  set  3  (topographic  heights)  to 
predict  values  of  additional  measured  gravity  anomalies  in  data 
set  2.  The  1983  Report  describes  the  work  of  three  sets  of 
investigators,  each  of  whom  used  some  variation  of  least- squares 
collocation.  Their  methods  and  results  are  summarized  briefly 
in  Table  3.2-1,  which  is  adapted  from  Table  2,  p.  92  of  the 
1983  Report.  The  standard  deviations  in  the  table  refer  to 
the  difference  between  the  predicted  anomalies  and  the  actual 

TABLE  3.2-1 

METHODS  FOR  COMPUTING  GRAVITY  ANOMALY 


METHOD 

AUTHORS 

STANDARD 
DEVIATION 
(mgal ) 

Leas  t- squares 
collocation  using 

Free-air  anomalies 

Merry 

Schwarz,  Lachapelle,  & 

17.8 

Mainville 

17.5 

Tscherning  &  Forsberg 

17.9 

Incomplete  Bouguer 

Merry 

7.4 

anomalies 

Schwarz,  Lachapelle,  & 

Mainville 

7.9 

Incomplete  Bouguer 
anomalies  plus 
trend  surface 

Merry 

6 . 6 

Refined  Bouguer 
anomalies 

Tscherning  &  Forsberg 

6.5 

Trend  surface  for 
incomplete  Bouguer 
anomalies 

Merry 

13.9 
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values  in  data  set  2.  These  results  show  that  least- squares 
collocation  (interpolation)  using  only  free-air  anomalies  is 
less  accurate  than  interpolation  based  on  Bouguer  anomalies. 


3.3  PREDICTION  OF  DEFLECTION  OF  THE  VERTICAL 

The  IAG  study  group  evaluated  algorithms  for  the  pre¬ 
diction  of  deflection  of  the  vertical  by  making  predictions 
based  on  free-air  gravity  in  data  sets  1  and  2,  supplemented 
by  the  gridded  topographic  heights  of  data  set  3.  The  results 
were  compared  to  the  astrogeodetic  deflections  of  data  set  4. 

Table  3.3-1  lists  six  algorithms.  The  first  four  are 
described  in  the  1983  Report  and  summarized  in  Ref.  46.  The 
last  four  are  described  and  systematically  compared  in  the 
1985  Report.  The  middle  two  overlap;  they  are  mentioned  in 
both  reports. 


TABLE  3.3-1 

METHODS  FOR  COMPUTING  DEFLECTION  OF  THE  VERTICAL 


METHOD 

AUTHORS 

Topographic/ Isostatic  Reduction 

Schwarz,  Lachapelle,  &  Mainville 
(1983  Report) 

Local  Collocation 

Hein  &  Landau  (1983  Report) 

Combined  Collocation  &  Integration 

Schwarz,  Lachapelle,  &  Mainville 
(1983  Report) ; 

Krynski  (1985  Report) 

Terrain  Effect  Integration 

Tscherning  &  Forsberg 

&  Collocation 

(1983,  1985  Reports) 

Rice  Rings  Integration 

Kearsley  (1985  Report) 

Fast  Fourier  Transform 

Sideris  (1985  Report) 

* 
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The  first  method,  Topographic/Isostatic  Reduction,  3 
the  only  one  in  which  the  gravity  data  are  not  used.  It  uses 
only  topography  and  the  low-degree  geopotential  model  GEM10B 
(Ref.  36).  An  Airy-Heiskanen  model  of  isostatic  compensation 
with  a  crustal  thickness  of  30  km  is  applied  and  the  effects 
of  the  topography,  isostatic  reduction  potential,  and  GEM10B 
are  combined  to  predict  deflection  of  the  vertical.  The  dif¬ 
ference  between  the  result  and  the  as trogeodetic  truth  data 
has  a  standard  deviation  ranging  from  2  to  2.5  arc  sec  (depend¬ 
ing  on  which  component  and  geographical  area  is  considered). 

The  next  method,  Hein  and  Landau's  pointwise  Local 
Collocation ,  is  applied  directly  to  the  free-air  gravity  anomaly 
data.  Values  within  a  radius  of  35  km  of  the  prediction  point 
are  used,  along  with  Reilly's  covariance  model  (e.g.,  Sec¬ 
tion  2.3.3)  fitted  to  the  data.  The  standard  deviations  rela¬ 
tive  to  the  truth  data  are  slightly  less  than  2  arc  sec. 

The  last  four  methods,  compared  in  detail  in  the  1985 
Report,  are  summarized  in  Table  3.3-2,  which  includes  a  brief 
discussion  of  each  method  and  an  indication  of  whether  deflec¬ 
tion  values  are  computed  separately  for  each  point  or  are  com¬ 
puted  simultaneously  (i.e.  ,  from  the  same  subset  of  the  gravity 
data)  for  more  than  one  point.  In  addition,  the  format  and  spa¬ 
tial  extent  of  the  input  gravity  data  and  the  treatment  of  local 
topographic  effects  and  global  effects  (from  long-wavelength 
geopotential  models)  are  indicated. 

Combined  Collocation  and  Integration  is  a  pointwise 
method  which  uses  (mean)  gravity  anomalies  within  a  3  deg  radius 
of  the  computation  point.  Least-squares  collocation  would  be 
inefficient  over  so  large  an  area  (too  much  data),  so  it  is 
applied  only  to  data  within  an  inner  zone  of  0.65  deg  radius, 
after  removal  of  the  long-wavelength  contribution  from  GEM10B . 
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The  Tscherning  and  Rapp  covariance  model  (Ref.  23;  Section  2.3.3), 
with  appropriate  parameters,  is  used  in  the  collocation.  Then, 
in  the  outer  annulus,  the  Vening  Meinesz  integral  is  applied 
to  the  difference  between  the  free- air  anomaly  data  and  the 
values  extrapolated  from  the  inner  zone  by  collocation.  The 
resulting  correction  is  applied  to  the  deflection  of  the  verti¬ 
cal  predicted  by  collocation,  and  the  effect  of  GEM10B  is  re¬ 
stored.  Finally,  an  additional  correction  is  applied  from  an 
approximate  Residual  Terrain  Model  (see  Ref.  34  and  the  follow¬ 
ing  paragraph) . 

Terrain  Effect  Integration  and  Collocation  is  also 
based  on  least-squares  collocation,  but  it  is  applied  over  a 
larger  area,  and  deflections  of  the  vertical  within  a  1  by  1  deg 
block  are  computed  simultaneously  (i.e.,  from  the  same  set  of 
gravity  anomalies)  rather  than  pointwise.  This  would  require 
excessive  computation  if  applied  to  the  raw  data;  however, 
prior  to  collocation,  the  effect  of  a  Residual  Terrain  Model 
(RTM)  is  removed  from  the  free-air  anomaly.  That  is,  the  ef¬ 
fect  of  topography  relative  to  a  smooth  reference  surface  is 
subtracted,  resulting  in  a  smooth  residual  field.  Also  removed 
is  the  effect  of  OSU79  (Ref.  37),  a  detailed  global  geopotential 
model  extending  to  spherical  harmonic  degree  and  order  180. 
Because  the  residual  field  is  smooth,  only  one  data  value  from 
each  3  by  3  arc  min  block  is  used  and  collocation  becomes  prac¬ 
tical.  The  covariance  function  used  in  the  collocation  is 
based  on  the  Tscherning  and  Rapp  covariance  model,  with  param¬ 
eters  suitable  to  the  residual  field.  After  collocation,  the 
RTM  and  OSU79  effects  are  restored  to  the  predicted  deflections. 

In  Rice  Rings  Integration  the  Vening  Meinesz  integral 
is  evaluated  by  means  of  a  computer  adaptation  of  Rice's  (Ref.  5) 
system  of  concentric  rings  divided  into  compartments  (see  Ap¬ 
pendix  A.l).  The  deflection  of  the  vertical  at  the  center  point 
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of  the  template  of  rings  is  computed  directly  from  the  free- 
air  anomaly  at  the  center  of  each  compartment.  In  turn,  the 
free-air  anomaly  at  the  center  of  a  compartment  is  estimated 
by  computing  the  incomplete  Bouguer  anomaly  at  nearby  points , 
applying  linear  interpolation,  then  using  the  height  at  the 
central  point  of  the  compartment  to  convert  back  from  Bouguer 
to  free-air  anomaly.  (This  is  the  only  way  in  which  the  gridded 
height  data  are  utilized.)  Finally,  the  deflection  predicted 
from  GEM10B  is  added  to  account  for  remote  zone  effects. 

The  Fast  Fourier  Transform  method  is  an  application 
of  the  FFT  approach  to  the  Vening  Meinesz  integral,  as  described 
in  Appendix  A. 4.  This  method  requires  gridded  data;  5  by  5  arc  min 
means  within  a  6  by  6  deg  block  are  used.  The  large  block  (the 
deflections  are  required  only  within  a  1  by  2  deg  area)  minimizes 
edge  effects.  The  complete  Bouguer  correction  (computed  via 
the  FFT  terrain  correction  algorithm  described  in  Appendix  B) 
and  the  effects  of  GEM10B  are  removed  from  the  input  free-air 
anomalies  and  restored  to  the  predicted  deflections.  The  FFT 
method  results  in  the  simultaneous  computation  of  deflections 
on  a  regular  grid,  which  must  be  interpolated  to  the  points  of 
interest . 


The  relative  performance  of  these  four  methods  is 
indicated  by  Table  3.3-3  (abstracted  from  Table  4.5,  p.  92, 

1985  Report).  The  statistics  in  the  second  column  from  the 
left  indicate  the  magnitude  and  mean  of  the  truth  data  (4  and  n  )  , 
while  the  last  four  columns  show  statistics  for  errors  (A£  and 
Aq)  in  the  prediction  of  each  algorithm.  The  statistics  include 
all  results  from  the  one  by  two  degree  area  between  32  and 
34  deg  North  and  106  and  107  deg  West.  Except  for  1.4  arc  sec 
in  the  case  of  An,  Rice's  Rings  Integration,  the  mean  errors 
are  less  than  1  arc  sec  in  absolute  value,  indicating  that 
there  is  little  systematic  bias.  The  standard  deviations  of 
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TABLE  3.3-3 

STATISTICS  FOR  COMPARISONS 


OBSERVED  MINUS  PREDICTED 

STATISTICS 

OBSERVED 
(Astrogeodetic ) 

COMBINED  COLLOCATION 
&  INTEGRATION 

TERRAIN  EFFECT 
INTEGRATION 
&  COLLOCATION 

RICE  RINGS 
INTEGRATION 

_ 

FAST  FOURIER 
TRANSFORM 

(arc  sec) 

4 

ISM 

A4  An 

A4 

40 

A4 

An 

A4 

An 

MEAN 

-2. A 

-A.  3 

-0.1  0.8 

-0.2 

— 

0.6 

-0.1 

1  .  A 

-0.1 

0.7 

STANDARD 

DEVIATION 

2.7 

6.2 

o.q  i.2 

0.9 

1.0 

0.9 

1.1 

1.0 

1.0 

RMS 

3.6 

7.5 

0.9  1 .  A 

0.9 

1.1 

_ 

0.9 

1.8 

1.0 

1.2 

the  errors,  including  both  deflection  components  and  all  four 
methods,  are  essentially  equal,  at  about  1  arc  sec.  Quite 
possibly,  this  represents  a  natural  limit,  imposed  by  the  spac¬ 
ing  and  quality  of  the  data,  which  has  been  reached  by  all 
four  algorithms. 

This  section  concludes  with  a  comment  on  relative 
efficiency  and  required  computer  time.  While  the  authors  of 
the  1985  Report  attempted  a  quantitative  comparison,  the  result 
is  not  straightforward.  This  is  because  the  data  preprocessing 
requirements  differ  among  the  algorithms,  and  the  relative  ef¬ 
ficiency  varies  with  the  number  of  computation  points.  The 
methods  which  involve  simultaneous  computations  (Terrain  Effect 
Integration  And  Collocation;  Fast  Fourier  Transform)  incur 
comparatively  less  computational  cost  for  generating  additional 
results  than  the  other  approaches.  Also,  the  terrain  correction 
computations  are  time  consuming  and  tend  to  dilute  other  dif¬ 
ferences.  However,  in  general,  the  Fast  Fourier  Transform 
method  tends  to  be  fastest  in  terms  of  computation  time  per 
estimated  deflection  point  and  the  other  three  algorithms  are 
roughly  comparable  to  one  another  in  speed. 
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4.  OVERVIEW  OF  TERRAIN  EFFECTS  STUDIES 

BASED  ON  WHITE  SANDS  DATA 


This  chapter  summarizes  TASC ' s  investigation  of  the 
significance  of  terrain  effects  in  the  computation  of  gravity 
disturbance  and  deflections  of  the  vertical,  complementing  the 
work  of  the  IAG  Study  Group  (Refs.  20,  21).  A  more  detailed 
description  of  this  study  is  found  in  Appendix  C. 

4.1  OVERVIEW 

Terrain  effects  need  to  be  included  in  computations 
that  predict  various  gravity  field  quantities  from  observed 
gravity  anomaly  data,  paticularly  in  regions  of  rugged  terrain 
like  the  White  Sands  study  area.  TASC ' s  investigation  of  the 
quantitative  significance  of  these  effects  is  outlined  in 
Fig.  4.1-1,  where  reference  is  made  to  the  sections  of  Appen¬ 
dix  C  in  which  detailed  technical  discussions  are  found,  and 
to  the  figures  in  this  appendix  which  present  specific  results 

Our  investigation  includes  two  distinct  studies.  One 
of  these,  shown  in  the  upper  half  of  Fig.  4.1-1,  makes  use  of 
spectral  techniques  to  examine  and  compare  observed  free-air 
gravity  anomalies  and  computed  topographic  effects  for  two 
types  of  terrain  --  mountain  and  plateau  --  and  for  two  pro¬ 
file  orientations  --  east-west  and  north-south.  The  methodol¬ 
ogy  and  conclusions  of  this  study  are  summarized  below  (Sec¬ 
tion  4.3)  and  presented  in  detail  in  Sections  C.2,  C.3,  and 
C.4  of  Appendix  C. 
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Figure  4.1-1  Overview  of  Terrain  Effects  Analyses 

The  second  study,  outlined  in  the  lower  half  of 
Fig.  4.1-1,  examines  how  well  deflections  of  the  vertical  can 
be  predicted  from  local  topographic  data  supplemented  by  a 
long-wavelength  component  derived  from  a  global  geopotential 
model.  This  study  is  also  summarized  below  (Section  4.4)  and 
is  documented  in  detail  in  Section  C.5  of  Appendix  C. 

The  results  and  conclusions  for  each  of  the  studies 
are  summarized  in  Section  4.5.  More  detail  will  be  found  in 
Appendix  C  (Sections  C.4.1,  C.4.2,  and  C.5. 2). 
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4.2  DATA  USED  IN  THESE  INVESTIGATIONS 


The  data  used  in  these  investigations  are  summarized 
in  the  left-hand  column  of  Fig.  4.1-1: 


•  A  set  of  3855  point  free-air  gravity 
anomaly  measurements,  covering  the  re¬ 
gion  from  32  to  35  deg  North  and  105  to 
108  deg  West  (White  Sands  study  area) 

•  Gridded  topographic  data  covering  the 
same  area,  sampled  at  a  density  of  120 
samples  per  degree 

•  A  set  of  geopotential  coefficients 
(Ref.  37)  used  to  characterize  long- 
wavelength  features  of  the  gravity  field 

•  A  set  of  441  deflections  of  the  vertical 
within  the  study  area,  observed  by  classi¬ 
cal  astrogeodetic  techniques,  used  as 
truth  data  to  evaluate  the  quality  of 
prediction  methods. 


The  free-air  gravity  anomaly  and  topographic  data  sets  are 
described  more  fully  in  Section  C.2  of  Appendix  C;  Fig.  C.2-1 
maps  the  locations  of  the  free-air  gravity  anomaly  data  points, 
while  Figs.  C.2-2  and  C.2-3  are  shaded  relief  maps,  in  color, 
of  the  gravity  anomalies  and  terrain  heights,  respectively. 
The  observed  deflection  data  are  described  in  Section  C.5  of 
Appendix  C  and  in  Ref.  20.  They  are  shown  graphically  in 
Fig.  C.5-1. 


4.3  SPECTRAL  ANALYSIS  OF  GRAVITY  DATA  AND 
TOPOGRAPHIC  EFFECT 

There  are  three  principal  steps  in  the  spectral  analy¬ 
sis  of  free-air  gravity  anomaly  data  and  computed  topographic 
effects : 
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•  Gridding  the  gravity  anomaly  data 

•  Selecting  profiles  and  computing  topo¬ 
graphic  effects  for  each  profile 

•  Carrying  out  spectral  analysis  to 
estimate : 

power  spectral  density  (PSD)  of 
the  gravity  anomalies  for  selected 
profiles 

power  spectral  density  (PSD)  of  the 
topographic  effect  for  selected 
profiles 

coherence  between  gravity  anomalies 
and  topographic  effect  for  selected 
subsets  of  the  profiles. 

Gridding  the  data  -  Since  the  free-air  gravity  anoma¬ 
lies  are  measured  at  points  unevenly  distributed  over  the  study 
area,  it  is  necessary  to  grid  the  data  before  the  spectral 
estimation  algorithms  can  be  applied.  The  gridding  algorithm 
is  described  in  Section  C.2  of  Appendix  C.  The  resulting  121  by 
121  grid  includes  40  data  points  per  degree.  Note  that  the 
topographic  data  are  already  supplied  in  regularly  gridded 
form.  Shaded  relief  maps,  in  color,  have  been  prepared  from 
the  gridded  anomaly  and  topography  data  (Figs.  C.2-2  and  C.2-3, 
Appendix  C. ) 

Selecting  profiles  -  Six  profiles  of  free-air  gravity 
anomaly  data  are  chosen  from  either  rows  or  columns  of  the 
gravity  anomaly  grid  described  above,  corresponding  to  east- 
west  and  north-south  orientations.  The  locations  of  the  three 
north-south  and  three  east-west  profiles,  selected  to  include 
mountain  and  plateau  terrain,  are  shown  in  Fig.  C.3-2  of 
Appendix  C. 
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Computing  topographic  effect  -  Topographic  effects 
are  estimated  by  computing  the  vertical  component  of  the  grav¬ 
ity  disturbance,  treating  terrain-height  data  as  a  grid  of 
vertically  oriented  rectangular  prisms  extending  upward  from  a 
sea-level  datum.  This  algorithm  is  described  in  Section  C.3 
of  Appendix  C.  The  computed  disturbances  are  related  to  the 
complete  Bouguer  correction  (Ref.  1). 

With  free-air  gravity  anomaly  and  topographic  effect 
computed  for  each  point  on  a  given  profile,  it  is  possible  to 
define  a  residual  anomaly  as  the  difference  between  the  two. 

An  example  is  given  in  Fig.  4.3-1  for  one  of  the  east-west 
profiles,  designated  as  EW-3.  This  figure  shows  the  variation 
of  terrain  elevation  along  the  profile,  the  observed  free-air 
gravity  anomalies,  the  computed  topographic  effect,  and  the 
residual  anomaly.  It  is  evident  from  the  figure  that,  for 
this  profile,  the  observed  free-air  gravity  anomalies  are  well 
correlated  with  topography,  but  that  the  residual  anomaly  has 
nearly  as  much  variability  as  the  observed  anomaly.  Similar 
figures  for  all  six  profiles  are  found  in  Appendix  C. 

Estimating  PSD  and  coherence  -  The  purpose  of  the 
power  spectral  density  computations,  described  in  detail  in 
Section  C.4  of  Appendix  C,  is  to  investigate  differences  be¬ 
tween  types  of  terrain  (mountainous  and  plateau)  and  to  examine 
the  significance  of  directional  anisotropy  (east-west  compared 
with  north-south).  In  order  to  make  spectral  comparisons,  the 
data  segments  are  pooled  into  four  groups: 

•  Plateau,  oriented  east-west 

•  Plateau,  oriented  north-south 

•  Mountain,  oriented  east-west 

•  Mountain,  oriented  north-south 


4-5 


THE  ANALYTIC  SCIENCES  CORPORATION 
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Figure  4.3-1  Profile  EW-3 


Power  spectra  for  each  of  the  four  groups  are  computed  by  state- 
space  modeling  (Ref.  44),  as  described  in  Appendix  C.  In  addi¬ 
tion,  the  relation  between  the  gravity  anomaly  and  topographic 
effect  data  series  is  summarized  in  a  coherence  plot,  which 
shows  the  degree  of  correlation  as  a  function  of  spatial  fre¬ 
quency.  PSD  and  coherence  plots  for  all  four  groups  are  pre¬ 
sented  and  discussed  in  Section  C.4  of  Appendix  C ;  one  example 
is  shown  here  as  Fig.  4.3-2. 
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In  this  example  of  mountainous  terrain  (east-west 
orientation),  the  free-air  anomaly  and  topographic  effect  have 
similar  spectra  (upper  plot  of  Fig.  4.3-2).  The  spectra  are 
nearly  flat  down  to  wavelengths  of  about  45  km.  The  coherence 
is  significant  over  a  broad  range  of  wavelengths,  showing  a 
peak  value  of  0.8  at  a  wavelength  of  about  40  km.  These  plots 
are  typical  of  east-west  profiles  (but  not  north- south ) ,  as 
discussed  in  Section  C.4.1  of  Appendix  C.  A  brief  summary  and 
interpretation  of  the  results  of  profile  spectral  analysis  is 
provided  in  Section  4.5. 

4.4  PREDICTING  DEFLECTIONS  OF  THE  VERTICAL  FROM 
LOCAL  TOPOGRAPHY 

There  have  been  a  number  of  investigations  (Refs.  20 
and  21,  for  example)  of  the  contribution  of  terrain  or  topo¬ 
graphic  effects  in  the  computation  of  deflections  of  the  verti¬ 
cal  from  gravity  anomaly  data.  This  phase  of  TASC ' s  study, 
documented  in  Section  C.5  of  Appendix  C,  examines  the  signifi¬ 
cance  of  topographic  effects  in  mountainous  areas  from  a  com¬ 
plementary  point  of  view.  Deflections  of  the  vertical  are 
computed  directly  for  441  locations  within  the  study  area  at 
which  observed  astrogeodetic  deflections  are  available,  using: 

•  Terrain  height  data  to  supply  local  or 
short-wavelength  information 

•  A  global  geopotential  model  to  supply 
long-wavelength  information. 

Gravity  anomaly  data  are  not  used.  The  deflections  predicted 
by  this  simple  method  are  then  compared  with  the  astrogeodetic 
truth  data.  The  results  are  discussed  in  Section  C.5.1  of 
Appendix  C . 
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A  graphical  overview  of  Che  relative  quality  of  this 
simple  prediction  method  is  given  in  Fig.  4.4-1,  which  shows 
the  observed  (left-hand  map)  and  predicted  (right-hand  map) 
deflections  in  vector-arrow  form.  The  length  and  direction  of 
an  arrow  indicate  magnitude  and  direction  of  the  deflection  of 
the  vertical;  the  tail  of  the  arrow  is  at  the  location  of  the 
deflection  station.  It  is  readily  observed  that  the  predicted 
deflections  correlate  well  with  the  direction  of  the  observed 
deflections,  but  the  predicted  magnitudes  tend  to  be  smaller. 

Another  view  of  the  correlation  between  observed  and 
predicted  deflections  of  the  vertical  is  given  by  scatter  dia¬ 
grams  like  the  one  shown,  as  an  example,  in  Fig.  4.4-2.  In 
this  diagram  the  predicted  value  of  the  east-west  component  of 
the  deflection  of  the  vertical  is  shown  as  the  x-component, 
and  the  corresponding  true  value  is  shown  as  the  y-component. 
Two  conclusions  may  be  drawn  from  Fig.  4.4-2: 

•  The  tight  clustering  of  the  points  ex¬ 
presses  a  high  coefficient  of  correla¬ 
tion  (approximately  equal  to  0.9) 

•  The  slope  of  the  least- squares  regression 
line  differs  from  one,  reflecting  a  system¬ 
atic  under-estimation  of  the  east-west 
component  of  the  deflection. 

Further  discussion  is  found  in  Appendix  C. 

4.5  SUMMARY  OF  RESULTS  AND  CONCLUSIONS 

For  the  first  investigation  reported  in  this  chapter, 
the  spectral  analysis  of  gravity  data  and  topographic  effect, 
the  results  are  presented  and  interpreted  in  Sections  C.4.1 
and  C.4.2  of  Appendix  C.  An  example  of  the  graphical  output 
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Figure  4.4-1  Observed  and  Predicted  Deflections  of  the  Vertical 
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Figure  4.4-2  Scatter  Diagram  of  East  Vertical 
Deflection  Component 


has  been  shown  in  Section  4.3  (Fig.  4.3-2).  The  principal 
conclusions  are: 


•  Anisotropy  (associated  with  the  regular 
north-south  strike  of  geologic  and  ter¬ 
rain  features  in  the  White  Sands  study 
area)  is  indicated  by  significant  dif¬ 
ferences  between  north-south  and  east- 
west  data  groups.  For  example,  there  is 
a  local  peak  in  the  PSDs  and  coherences 
in  the  east-west  groups,  but  not  in  the 
north-south  groups. 
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•  Mountain  and  plateau  areas  differ  in 
several  ways,  among  which  are: 

Coherence  between  topographic  ef¬ 
fects  and  free-air  gravity  anoma¬ 
lies  is  low  at  long  wavelengths 
for  plateau  areas,  but  high  at 
long  wavelengths  in  the  mountains 

In  plateau  areas  the  power  in  the 
topographic  effect  is  lower  than 
the  power  in  the  free-air  and  re¬ 
sidual  anomalies,  while  for  moun¬ 
tain  areas  the  residual  anomaly 
has  the  lowest  power  at  long 
wavelengths . 


All  of  these  differences  are  consistent  with  the  anisotropic 
and  nonstationary  nature  of  the  terrain  and  gravity  anomalies 
in  the  White  Sands  region. 

The  results  of  the  second  phase  of  TASC ' s  investiga¬ 
tion,  the  prediction  of  deflections  of  the  vertical  from  local 
topographic  data,  are  summarized  and  discussed  in  Section  C.5.1 
of  Appendix  C.  Examples  of  the  results,  in  graphical  form, 
have  been  shown  in  Section  4.4  (Figs.  4.4-1  and  4.4-2).  The 
major  conclusion  is  a  conf irmation  of  the  importance  of  account¬ 
ing  for  local  topographic  effects  in  the  computation  of  deflec¬ 
tions  of  the  vertical,  since  --  in  this  predominantly  mountainous 
region  --  the  topographic  effect  by  itself  accounts  for  most 
of  the  deflection  signal. 
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5.  NEW  DERIVATION  OF  THE  FFT  TERRAIN 

CORRECTION  ALGORITHM 


The  gravimetric  terrain  correction  is  important  in 
mountainous  areas,  but  is  very  time  consuming  to  compute  in 
its  exact  form,  since  it  is  a  nonlinear  function  of  the  ter¬ 
rain  elevation  h(x,y).  Recently,  an  approximate  form  of  the 

terrain  correction,  in  terms  of  two-dimensional  linear  con- 

2 

volution  integrals  involving  h  and  h  ,  has  been  derived 
(Refs.  33  and  34).  The  convolutions  may  be  performed  effi¬ 
ciently  using  gridded  terrain  data  and  the  Fast  Fourier  Trans¬ 
form  (FFT),  by  virtue  of  the  convolution  theorem  (e.g.,  Appen¬ 
dix  A. 4).  This  approach  has  been  successfully  applied  to  com¬ 
pute  the  complete  Bouguer  anomaly  in  the  White  Sands  area 
(Ref.  21),  for  example.  However,  the  published  derivations 
include  ad  hoc  steps  that,  apparently,  are  required  to  ensure 
good  numerical  behavior. 

This  chapter,  in  contrast,  describes  the  essentials 
of  a  new  derivation  (Ref.  35)  of  the  approximate  terrain  cor¬ 
rection  that  avoids  such  ad  hoc  steps.  It  shows  a  more  rigorous 
and  direct  route  to  the  same  final  result  and  clarifies  the 
physical  meaning  and  optimum  choice  for  a  length  parameter 
that  was  introduced  arbitrarily  in  Refs.  33  and  34.  The  pres¬ 
entation  here  is  brief  and  descriptive;  further  details  are 
given  in  Appendix  B. 

The  terrain  correction  is  to  be  added  to  the  gravity 
anomaly  following  subtraction  of  the  simple  Bouguer  (slab) 
correction.  At  a  point  P,  it  corresponds  to  the  effect  of  the 
true  topography  minus  a  slab  of  height  hp,  where  hp  =  h  (Xp,Yp) 
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is  the  elevation  of  P.  This  is  shown  schematically  in  Fig.  5-1. 
Above  the  computation  point  P  there  are  hills  of  excess  mass 
and  below  P  there  are  valleys  of  (effectively)  negative  mass. 
Both  have  the  effect  of  decreasing  the  observed  gravity  anomaly; 
the  terrain  correction  is  always  positive  *  Assuming  a  constant 
density  p  for  the  topography  and  denoting  the  gravitational 
constant  by  G,  the  exact  (flat-earth)  terrain  correction  is 


h(x,y) 


c(xp,yp)  =  Gp  If  dxdy 


3 if  dX<iy  •/ 


(z  -  hp) 


dz 


(xp-x)2  +  (yp-y)2  +  (hp-z)2 
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(5-1) 
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COMPUTATION  POINT 
2  =  h  (Xp.  Yp)  =  hp 


TERRAIN  SURFACE: 
Z  =  h  (X,  Y) 
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Starting  from  Eq.  5-1,  previous  derivations  (Refs.  33 
and  34)  of  the  FFT  terrain  correction  algorithm  have  consisted  of 

1.  An  approximation  leading  to  convolutions  of  h 
2 

and  h  with  the  singular  kernel 


2  2 
(Xz  +  Y  ) 


-3/2 


which  is  not  integrable  over  the  x-y  plane  and 
whose  Fourier  transform  is  nonexistent 

Replacement  of  the  singular  kernel  by 


(X*  +  Y  +  a  ) 


-3/2 


which  has  the  Fourier  transform 


—  exp  (-2naVu^+v^) 
3 


(following  the  convention  of  Appendix  A. 4  and 
Ref.  32). 


(In  Ref.  33  an  alternate  but  equally  arbitrary  version  of  step  2 
is  also  applied.) 

In  fact,  however,  convolution  of  h(x,y)  and  its  square 

with  the  modified  kernel  is  not  only  more  tractable  numerically, 

2  2  -3/2 

but  also  more  accurate  than  convolution  with  (X  +Y  )  ,  for 

a  range  of  nonzero  values  of  a.  This  result  follows  directly 

from  the  exact  formula  (Eq.  5-1)  by  means  of  a  single  approxima- 

2  2 
tion  in  which  (hp-z)  is  replaced  by  a  constant  a  , 


h(x,y) 


c(xp,yp)  =  Gp 


(Xp-x 


(z  -  hp)  dz 
,2  ~  7T~  213/2 

)  +  (yp-y)  ♦  a  | 
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Starting  from  Eg .  5-2,  integration  over  z  and  expansion  of  the 

2  2 

numerator  lead  directly  to  convolution  of  h  and  h  with  (X  + 

Y2+a2)'3//2  (for  details,  see  Appendix  B).  This  yields  the 

same  result  as  steps  1  and  2  above,  but  now  it  is  clear  that 
2  o 

ar  represents  an  average  value  of  (hp-z)  and  that  its  optimal 

value  always  exceeds  zero.  The  two  steps  of  approximation 
applied  in  Refs.  33  and  34  are  unnecessary  and  partially  can¬ 
cel  one  another. 

As  a  demonstration  that  the  best  value  of  a  is  indeed 
nonzero,  Fig.  5-2  shows  a  plot  of  rms  error  in  the  terrain 
correction  vs  the  value  of  a  for  a  sample  of  100  points  in  a 
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rugged  portion  of  the  White  Sands  area  33.1  to  33.5  deg  North 
and  105.5  to  105.9  deg  West.  The  error  was  computed  by  sub¬ 
tracting  the  approximate  terrain  correction  (Eq.  5-2),  com¬ 
puted  in  terms  of  convolutions,  from  the  result  of  evaluating 
the  exact  formula  (Eq.  5-1)  by  numerical  (bicubic  spline)  inte¬ 
gration.  It  is  clear  that  the  best  value  of  a  for  the  study 
case  is  about  250  m,  and  that  any  choice  of  a  up  to  an  upper- 
bound  of  about  375  m  will  yield  a  lower  rms  error  than  a  =  0. 
The  rms  errors  shown  are  quite  small,  but  the  maximum  absolute 
errors  are  larger;  0.84  mgal  for  a  =  0  vs.  0.28  mgal  for 
a  =  250  m.  (The  appropriate  choice  of  a  leads  to  a  factor- 
of-3  improvement.) 

In  summary,  this  chapter  and  the  more  detailed  deri¬ 
vation  and  discussion  in  Appendix  B  show  that  the  FFT  terrain 
correction  algorithm  is,  in  fact,  more  rigorous  than  indicated 
by  its  proponents  in  Refs.  33  and  34,  and  that  the  optimal 
value  of  the  parameter  a  always  exceeds  zero. 
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6.  NUMERICAL  EXPERIMENTS  WITH  GEOFAST  ALGORITHM 

6.1  SUMMARY  OF  APPROACH 

The  GEOFAST  Algorithm  is  described  in  Appendix  A. 3. 

This  technique  computes  minimum-variance  estimates  of  geodetic 
quantities  (equivalent  to  least-squares  collocation)  in  a  highly 
efficient  manner.  The  algorithm  accepts  two-dimensional  gridded 
data,  and  takes  advantage  of  the  structure  of  stationary  statis¬ 
tical  gravity  models.  High  computational  efficiency  is  achieved 
by  frequency-domain  techniques  based  on  the  Fast  Fourier  Trans¬ 
form  (FFT).  This  approach  yields  a  computational  load  propor¬ 
tional  to  N  log  N,  where  N  is  the  number  of  data  points,  as 
3 

compared  to  N  for  the  matrix  inversion  used  in  standard  collo¬ 
cation.  Appendix  A. 3  describes  how  the  covariance  matrices  of 
the  collocation  method  are  transformed  into  the  frequency  domain 
by  a  two-dimensional  windowed  FFT.  The  resulting  frequency- 
domain  covariance  matrices  are  approximately  block-banded,  in 
the  sense  that  the  out-of-band  elements  are  below  a  threshold,  e. 

There  is  a  tradeoff  between  the  number  of  bands  retained,  mD  , 

o 

and  the  accuracy  of  the  approximation,  e.  The  computer  load 

2 

for  the  GEOFAST  solution  is  asymptotical ly  proportional  to  m_ 

B 

as  illustrated  in  Fig.  6.1-1  for  a  32  by  32  grid  (N  =  1024) 
using  the  VAX  11/780  computer.  For  the  case  shown,  a  standard 
collocation  solution  would  require  more  than  two  hours  (vs 
15  min  for  GEOFAST  with  m^  =  8). 

Two  sets  of  numerical  experiments  were  conducted  with 
the  GEOFAST  algorithm.  The  first  set  used  synthetic  marine 
gravity  data  in  an  ocean  trench  area.  This  data  set  was  gen¬ 
erated  using  the  Rapp  1981  spherical  harmonic  expansion  to 
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Figure  6.1-1  GEOFAST  Computer  Time  for  N  =  1024 

degree  ai»u  order  180  (Ref.  38).  Such  a  model  resolves  wave¬ 
lengths  down  to  about  200  km.  The  second  set  of  experiments 
used  a  higher  frequency  local  model  developed  by  TASC  and  NSWC 
for  gravity  gradiometer  evaluation  (Ref.  39).  This  data  set 
is  based  on  a  mass  dipole  realization  and  contains  wavelengths 
as  short  as  5  km. 

The  experiments  consisted  of  using  the  GEOFAST  algo¬ 
rithm  to  compute  deflections  of  the  vertical  from  gravity  anoma¬ 
lies  for  each  of  the  two  data  sets.  Parameters  which  were 
explored  include  the  data  extent,  grid  spacing,  GEOFAST  band¬ 
width,  covariance  models,  and  measurement  noise  level.  A  de¬ 
scription  of  the  two  data  sets  and  the  results  of  the  experi¬ 
ments  are  given  in  the  following  sections. 
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6.2  SYNTHETIC  MARINE  GRAVITY  DATA  SET 

This  data  set  was  generated  from  the  Rapp  1981  world¬ 
wide  spherical  harmonic  expansion  to  degree  and  order  180 
(Ref.  38).  The  test  area  selected  is  a  36  deg  square  in  the 
Western  Pacific  which  includes  a  significant  ocean  trench 
(Fig.  6.2-1).  This  model  expansion  is  capable  of  resolving 
only  wavelengths  longer  than  about  200  km.  Both  gravity  anoma 
lies  and  deflections  of  the  vertical  are  generated  for  grid 
spacings  of  0.5  deg  and  0.25  deg.  The  smaller  spacing  repre¬ 
sents  "oversampling"  but  is  included  to  study  algorithm  depend 
ence  on  the  number  of  datapoints.  The  data  sets  obtained  con¬ 
stitute  the  "truth  data"  for  the  experiments  described  below. 


LONGITUDE  (E) 

Synthetic  Marine  Gravity  Anomalies  (Rapp 
Model  to  Degree  180  Contoured  at  20  mgal) 


Figure  6.2-1 
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To  study  the  effect  of  high-order  reference  fields  on 
gravity  estimation  techniques,  the  Rapp  model  is  also  used  to 
generate  anomalies  and  deflections  for  spherical  harmonic  de¬ 
grees  and  orders  91  through  180.  These  data  (Fig.  6.2-2)  cor¬ 
respond  to  subtracting  a  degree  90  reference  field  from  the 
data  prior  to  GEOFAST  processing.  The  resulting  data  set  con¬ 
tains  only  the  high-frequency  portion  of  the  field  (in  this 
case,  wavelengths  between  200  and  400  km).  The  removal  of  the 
reference  field  increases  the  accuracy  of  the  estimates.  It 
is  compensated  for  by  restoring  the  reference  field  to  the 
estimated  quantities. 


A -40257 


Figure  6.2-2  High-Frequency  Synthetic  Marine  Gravity 
Anomalies  (Rapp  Model  Degrees  91-180 
Contoured  at  20  mgal) 
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6.3  RESULTS  WITH  SYNTHETIC  MARINE  DATA 
6.3.1  Description  of  Test  Cases 

The  first  three  test  cases  use  truth  data  from  the 
high-frequency  Rapp  model  (degrees  91  through  180).  The  ex¬ 
tent  of  the  data  is  reduced  to  the  northwest  15  deg  square  in 
Fig.  6.2-2  (latitude  17  deg  to  32  deg  North,  and  longitude 
135  deg  to  150  deg  East).  Since  the  trench  runs  north-south, 
east  deflections  of  the  vertical  were  estimated  from  gravity 
anomalies.  The  truth  data  (both  anomalies  and  deflections) 
were  gridded  at  0.5  deg  (Figs.  6.3-1  and  6.3-2)  and  at  0.25  deg 
(Fig.  6.3-3).  The  resulting  data  sets  contain  900  points  (30 
by  30)  and  3600  points  (60  by  60),  respectively.  The  gravity 
anomaly  in  this  region  varies  from  -109  mgal  to  137  mgal ,  while 
the  east  deflection  varies  from  -19  arc  sec  to  25  arc  sec. 

The  r ms  value  for  the  east  deflection  is  6  arc  sec. 


150.0  '  17.0 

Figure  6.3-1  Synthetic  Gravity  Anomaly  Data  (Rapp  Model 
Degrees  91-180  at  0.5  deg  Spacing) 


6-5 


THE  ANALYTIC  SCIENCES  CORPORATION 


THE  ANALYTIC  SCIENCES  CORPORATION 


The  covariance  model  for  the  GEOFAST  estimation  was  the 
Attenuated  White  Noise  (AWN)  model  developed  at  TASC  (Ref.  24). 
This  model  has  several  desirable  properties,  especially  the 
availability  of  closed-form  expressions  for  all  related  grav¬ 
ity  quantities,  and  the  ability  to  closely  fit  empirical  co- 
variance  data  (see  Appendix  A. 2. 3).  For  this  test,  a  single¬ 
shell  AWN  model  was  chosen  with  a  correlation  distance  of 
0.5  deg.  It  is  not  necessary  to  specify  the  gravity  model 
field  strength  (variance)  since  this  scale  factor  cancels  out 
in  the  estimation  gains.  For  the  initial  test  cases  the  ef¬ 
fect  of  measurement  noise  (uncertainty  in  the  anomaly  data) 
was  ignored  and  the  exact  Rapp  model  data  were  used  for  esti¬ 
mation.  Degradation  in  the  estimates  due  to  measurement  noise 
is  analyzed  separately  by  a  perturbation  approach  and  is  dis¬ 
cussed  in  Section  6.3.3. 

6.3.2  Results  for  Noise-Free  Estimation 

Table  6.3-1  displays  the  results  of  three  test  cases: 

•  0.25  deg  grid  with  mg  =  0 

•  0.25  deg  grid  with  mg  =  4 

•  0.50  deg  grid  with  mg  =  0  . 

Estimation  error  is  given  in  units  of  arc  sec  for  the  central 
12  deg  square  (within  the  15  deg  square  of  data).  This  re¬ 
striction  reduces  the  influence  of  edge  error  due  to  data 
truncation,  which  is  common  to  all  estimation  methods.  The 
deflection  estimates  are  within  0.2  arc  sec  (standard  devi¬ 
ation)  for  all  cases  and  within  0.7  arc  sec  (absolute  maximum). 
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TABLE  6.3-1 

GEOFAST  ESTIMATION  ACCURACY 


RESOLUTION 

EFFICIENCY 

ESTIMATION  ERROR 

irk 

(arc  sec) 

GRID  (deg) 

POINTS 

“b 

TIME  (min)* 

MEAN 

SIGMA 

MAXIMUM 

0.25 

3600 

0 

4 

0.001 

0.15 

0.70 

0.25 

3600 

4 

30 

0.003 

0.09 

0.31 

0.50 

900 

0 

1 

0.003 

0.17 

0.59 

*Computer  time  for  VAX  11/780 
**Central  12  deg  square 

As  expected,  results  are  best  at  the  finest  grid  and 
largest  bandwidth  --  the  rms  and  maximum  errors  are  reduced  by 
a  factor  of  two  (0.09  arc  sec  and  0.31  arc  sec,  respectively) 
compared  to  the  coarse  grid  and  zero  bandwidth.  However,  the 
effect  of  reducing  the  grid  spacing  alone  (at  zero  bandwidth) 
is  only  a  modest  improvement  in  accuracy.  This  result  is  due 
to  the  oversampling  (0.25  deg  is  about  28  km)  when  the  data 
wavelengths  are  greater  than  200  km.  The  mean  error  is  small 
(less  than  0.003  arc  sec)  in  all  cases  because  of  removal  of 
the  high  order  reference  field  (Rapp  model  to  degree  90). 

Figure  6.3-4  gives  a  detailed  plot  of  the  estimation 
error  at  each  grid  point  for  two  of  the  cases.  In  general, 
the  error  is  largest  near  the  edges  or  where  the  field  is 
largest,  and  varies  from  -0.77  arc  sec  to  0.67  arc  sec  for  the 
best  case. 

The  efficiency  of  the  GEOFAST  algorithm  is  also  evident 
from  Table  6.3-1.  Deflections  accurate  to  0.2  arc  sec  rms  are 
obtained  in  1  min  of  VAX  11/780  computer  time.  The  table  also 
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shows  the  increase  in  computer  time  to  four  minutes  when  the 

grid  density  is  doubled.  Increasing  the  bandwidth  has  the 

largest  influence  on  computer  loading,  raising  computing  time 

for  3600  data  to  30  minutes  when  M„  =  4.  These  results  demon- 

o 

strate  the  tradeoff  available  through  choice  of  the  bandwidth 
to  achieve  a  desired  accuracy  level  commensurate  with  the  in¬ 
herent  data  quality.  In  the  standard  collocation  approach, 
the  user  is  forced  to  accept  (and  pay  for)  full  numerical  ac¬ 
curacy  which  substantially  exceeds  the  confidence  limits  set 
by  measurement,  gridding,  and  truncation  errors.  Note  that 
the  computer  time  required  for  standard  collocation  using  3600 
data  would  be  greater  than  GEOFAST  by  at  least  two  orders  of 
magnitude.  (An  exact  determination  of  this  is  not  practical 
as  collocation  runs  using  more  than  500  data  at  once  are  pro¬ 
hibitively  expensive.) 

Effect  of  Reference  Field  Removal  --  To  quantify  the 
benefits  of  removing  a  high  order  reference  field  from  the 
data  (as  was  done  in  the  three  cases  discussed  above),  the 
first  test  case  (0.25  deg  grid  and  mg  =  0)  was  repeated  using 
data  from  the  full  Rapp  model  (through  degree  and  order  180). 
The  results  are  shown  in  Table  6.3-2.  Note  that  removal  of 
the  reference  field  reduces  the  rms  field  strength  from  9.2  arc 
sec  to  6.1  arc  sec.  The  results  show  a  sizeable  mean  error 
(-1.5  arc  sec)  when  no  reference  field  is  removed.  This  ef¬ 
fect  results  from  the  lack  of  low-frequency  (long  wavelength) 
information  in  local  data  and  could  have  been  reduced  substan¬ 
tially  by  removing  a  low-order  reference  field  (e.g.,  to  degree 
and  order  12).  However,  when  no  reference  field  is  removed, 
the  degradation  in  standard  deviation  (by  a  factor  ot  5)  and 
maximum  error  (by  a  factor  of  nearly  10)  are  quite  dramatic 
and  justify  the  use  of  the  high-order  reference  field. 
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TABLE  6.3-2 

EFFECT  OF  HIGH-ORDER  REFERENCE  FIELD 


REFERENCE 

FIELD 

REMOVED 

REMAINING 
FIELD  rms 
(arc  sec) 

ESTIMATION  ERROR* 

(arc  sec) 

MEAN 

SIGMA 

MAXIMUM 

None 

9.16 

-1.520 

0.690 

2.960 

Degree  90 

6.14 

0.001 

0.120 

0.310 

^Central  5  degree  square 


For  any  geodetic  algorithm,  estimation  errors  due  to 
finite  extent  of  data  decrease  towards  the  center  of  the  esti¬ 
mation  region.  Thus,  the  estimation  errors  shown  in  the  second 
row  of  Table  6.3-2  (based  on  the  5  deg  square)  are  slightly 
smaller  than  the  results  in  the  first  row  of  Table  6.3-1  (based 
on  the  12  deg  central  square  of  the  same  GEOFAST  output). 

6.3.3  Effect  of  Noisy  Observations  on  Estimation  Results 


The  results  in  Table  6.3-1  and  Table  6.3-2  are  based 
on  perfect  (error- free)  anomaly  data,  and  thus  are  only  a  lower 
bound  on  attainable  accuracy.  The  effect  of  measurement  er¬ 
rors  on  these  results  was  assessed  analytically,  since  the 
data  errors  are  normally  a  small  fraction  of  the  field  strength 
(high  signal-to-noise  ratio).  The  analysis  is  based  on  the 
Wiener  smoother  formulation  for  continuous  data  in  an  infinite 
plane  (Refs.  40,  41),  except  that  the  gridding  (discrete  sam¬ 
pling)  is  accounted  for.  The  gravity  field  is  assumed  to  be 
band  limited  and  the  sample  spacing  to  satisfy  the  Nyquist 
criterion.  As  a  result,  there  is  no  aliasing  error  due  to 
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sampling.  These  conditions  are  satisfied  here  since  the  data 
in  the  Rapp  model  contain  no  spectral  energy  at  frequencies 
higher  than  the  Nyquist  frequency  (1/2L  where  L  is  the  grid 
spacing) . 


The  Wiener  transfer  function  (between  east  deflection 
and  anomaly) ,  corresponding  to  the  AWN  model  for  the  data  plus 
white  noise,  is  easily  obtained  in  the  frequency  domain.  To  a 
first  (linear)  approximation,  the  effect  of  the  noise  is  addi¬ 
tive  at  the  variance  level  and  can  be  calculated  separately  as 
a  "noise-induced"  error  in  the  estimates.  This  error  is  com¬ 
bined  with  the  noise- free  GEOFAST  estimation  error  (from 
Table  6.3-1)  to  give  the  total  estimation  error  for  noisy  data 
The  result  of  this  calculation  for  various  noise  levels  is 
shown  in  Table  6.3-3.  For  example,  for  5  mgal  anomaly  errors 
on  a  0.25  deg  grid,  the  effect  is  to  degrade  the  total  estima¬ 
tion  error  from  0.15  arc  sec  to  0.22  arc  sec  rms . 

TABLE  6.3-3 

EFFECT  OF  MEASUREMENT  ERROR 


GRID  (deg) 

MEASUREMENT  ERROR 
rms  (mgal) 

NOISE- INDUCED  ERROR 
rms  (arc  sec) 

TOTAL 

ESTIMATION  ERROR 
rms  (arc  sec) 

0.25 

0 

0.00 

0.15 

1 

0.03 

0.15 

2 

0.07 

0.  17 

5 

0.  16 

0.22 

0.50 

0 

0.00 

0.17 

1 

0.04 

0.17 

2 

0.08 

0.19 

5 

0.20 

0.26 
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6.4  SYNTHETIC  LOCAL  GRAVITY  DATA  SET 

This  data  set  was  generated  for  use  in  the  gravity 
gradiometer  testing  program  (Ref.  39).  A  mass  dipole  model 
was  constructed  with  prescribed  statistical  properties.  The 
statistical  gravity  model,  representative  of  the  gradiometer 
test  area  in  north  Texas,  is  specified  by  a  seven-shell  AWN 
gravity  potential  model  as  shown  in  Table  6.4-1.  The  result¬ 
ing  value  of  rms  vertical  disturbance  (equivalent  to  the  grav¬ 
ity  anomaly  for  a  flat-earth  approximation)  is  31.4  mgal . 

This  model  has  significant  spectral  power  at  wavelengths  of 
10  km  (Fig.  6.4-1) . 

Synthetic  data  were  generated  from  the  mass  dipole 
model  for  the  vertical  gravity  disturbance  (equivalent  to  grav¬ 
ity  anomaly)  and  the  horizontal  gravity  disturbances  (equiva¬ 
lent  to  deflections  of  the  vertical).  These  data  consist  of  a 
500  km  square  gridded  at  5  km  spacing  and  constitute  the  truth 
data  for  the  experiments  described  below.  As  with  the  syn¬ 
thetic  marine  gravity  data,  the  GEOFAST  algorithm  is  used  to 
estimate  the  east  gravity  disturbance  from  the  vertical  gravity 
disturbance . 


TABLE  6.4-1 

AWN  STATISTICAL  GRAVITY  MODEL 


SHELL  NO. 

DEPTH 

(km) 

RMS  POTENTIAL 
(mgal • km) 

1 

2.1 

2.3 

2 

5.0 

11.0 

3 

16.0 

72.0 

4 

52.0 

580.0 

5 

161.0 

2300.0 

6 

861.0 

7000.0 

7 

2150.0 

33000.0 
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Figure  6.4-1  Along-Track  Vertical  Disturbance  Spectra 


Figure  6.4-2  displays  a  surface  plot  of  a  150  km 
square  region  of  the  east  gravity  disturbance  data  at  5  km 
spacing  (900  points).  A  contour  map  of  the  same  data  is  pre¬ 
sented  in  Fig.  6.4-3.  There  is  a  pronounced  west-to-east 
trend  from  positive  to  negative  values  in  the  range  36  mgal  to 
-25  mgal.  The  average  value  of  the  east  disturbance  is 
1.03  mgal  and  the  rms  is  12.3  mgal.  (This  is  equivalent  to 
about  2.5  arc  sec  of  deflection.) 

6.5  RESULTS  WITH  SYNTHETIC  LOCAL  DATA 

One  test  case  was  computed  with  truth  data  from  the 
high-frequency  local  gravity  data  sets.  No  reference  field 
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was  removed  from  the  data.  The  extent  of  the  data  was  a  150  km 
square  gridded  at  5  km  for  a  total  of  900  points.  The  covari¬ 
ance  model  for  the  GEOFAST  estimation  was  a  single- shell  AWN 
model  with  a  correlation  distance  of  10  km.  No  measurement 
noise  was  included  since  this  effect  can  be  assessed  analyti¬ 
cally  as  in  Section  6.3.3.  The  GEOFAST  bandwidth  parameter 
was  chosen  to  be  m^  =  4 . 

Estimation  accuracy  results  for  this  case  are  dis¬ 
played  in  Table  6.5-1  for  successively  smaller  estimation  re¬ 
gions  within  the  dataset.  There  is  a  systematic  negative  bias 
of  about  7  mgal  which  represents  the  low-frequency  field  not 
observable  in  such  a  small  data  set  (150  km  is  less  than  2  deg). 
This  effect  would  be  removed  by  first  subtracting  the  mean  or 
by  a  sufficiently  high  order  reference  field,  as  demonstrated 
in  the  previous  section.  The  error  standard  deviation  decreases 
monotonically  toward  the  center  of  the  data  (from  almost  9  mgal 
to  less  than  2  mgal).  The  accuracy  near  the  center  is  1.7  mgal 
or  equivalently  0.35  arc  sec  of  deflection.  The  rms  east  dis¬ 
turbance  field  in  the  same  region  is  9.5  mgal  so  that  the  rms 


TABLE  6.5-1 

GEOFAST  ESTIMATION  ACCURACY 


*Side  length  of  central  square  within  150  km 
square  of  data. 
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within  1  or  2  correlation  distances  from  the  edge.  Removal  of 
a  high-order  reference  field  would  reduce  this  effect  by  reduc¬ 
ing  the  correlation  distance  inherent  in  the  data.  In  addi¬ 
tion,  some  of  the  sharp  central  valleys  in  Fig.  6.4-2  have  not 
been  resolved  precisely  enough.  For  this  data  set  a  finer 
grid  spacing  would  be  expected  to  provide  significant  improve¬ 
ment  here,  by  more  precisely  aligning  estimated  peaks  with  the 
truth  data. 


6.6  SUMMARY  OF  RESULTS 

Numerical  experiments  were  conducted  with  the  GEOFAST 
algorithm  using  synthetic  data  from  two  different  sources. 
The  first  data  set  was  the  Rapp  1981  spherical  harmonic  model 
evaluated  through  degree  and  order  180  for  a  marine  trench 
area  in  the  Pacific.  This  provides  a  moderately  high  fre¬ 
quency  model  in  an  active  ocean  gravity  area. 

Numerical  experiments  with  the  first  data  set  demon¬ 
strate  that  GEOFAST  is  able  to  estimate  deflections  of  the 
vertical  from  gravity  anomalies  to  less  than  0.2  arc  sec  rms 
error.  This  result  is  dependent  on  reducing  the  data  by  a 
high-order  reference  field  (in  this  case  the  Rapp  model  to 
degree  90).  With  only  a  mean  subtracted  from  the  data,  esti¬ 
mation  accuracy  is  degraded  to  the  order  of  1.0  arc  sec  rms  (a 
f ac tor-o f - f ive  degradation). 

The  effect  of  measurement  noise  in  anomaly  data  on 
deflection  estimation  accuracy  was  also  studied.  For  white 
noise  up  to  5  mgal  rms,  total  estimation  errors  were  less  than 
0.3  arc  sec  rms. 
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The  efficiency  of  the  GEOFAST  algorithm  was  demon¬ 
strated  by  quadrupling  the  number  of  datapoints ,  N,  from  900 
to  3600.  The  cost  per  calculation  increased  as  predicted  pro¬ 
portional  to  N  log  N.  In  contrast,  standard  collocation  ap- 

3 

proaches  have  a  computer  load  proportional  to  N  .  Another 
important  feature  of  the  GEOFAST  algorithm  that  is  also  illus¬ 
trated  is  the  design  tradeoff  between  computing  cost  and  esti¬ 
mation  accuracy  (through  a  window  bandwidth  parameter).  By  this 
technique,  a  desired  level  of  accuracy  can  be  achieved  for  the 
minimum  computer  load.  In  the  cases  reported  here,  relaxing 
the  accuracy  requirement  from  0.1  arc  sec  to  0.2  arc  sec  re¬ 
sults  in  a  reduction  of  computation  time  by  a  factor  of  8. 

The  second  source  of  synthetic  data  employed  in  these 
experiments  is  a  high-frequency  local  gravity  model  at  a  reso¬ 
lution  of  5  km  in  a  moderately  active  area  of  the  U.S.  The 
GEOFAST  algorithm  was  used  with  this  data  set  to  predict  hori¬ 
zontal  gravity  disturbance  from  vertical  gravity  disturbance 
(equivalent  to  deflection  prediction  from  anomalies).  Results  \ 

show  that  estimation  accuracy  better  than  2  mgal  rms  (0.5  arc 
sec  deflection)  is  attainable  away  from  the  edge  of  the  data. 

This  level  of  accuracy  requires  subtracting  out  a  high-order 
reference  field  to  eliminate  wavelengths  on  the  order  of  100  km 
from  the  data.  To  obtain  additional  accuracy,  this  gravity 
field  would  need  to  be  sampled  at  a  higher  rate  than  5  km. 

The  results  presented  in  this  chapter  demonstrate 
that  the  GEOFAST  technique  could  be  applied  to  current  prob¬ 
lems  in  reduction  of  gravity  gradiometer  or  satellite  altimeter 
data.  In  these  applications,  a  large  volume  of  gridded  data 
is  available  and  computer  loading  poses  the  major  bottleneck. 

Accuracy  and  resolution  requirements  inherent  in  high-quality 
map  production  dictate  a  need  to  process  the  data  without  com¬ 
pression  via  averaging  or  parametric  surface  fitting.  In  this 
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context,  the  computer  problem  must  be  addressed  directly,  and 
the  GEOFAST  method  is  an  example  of  such  an  approach.  In  con¬ 
trast,  the  collocation  approaches  described  in  Refs.  20  and  21 
and  summarized  in  Chapter  3  were  unable  to  utilize  all  of  the 
available  data  observations.  These  approaches  required  aver¬ 
aging  or  subsampling  the  observations  (see  Table  3.3-2)  before 
estimation . 
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7. 


SUMMARY  AND  CONCLUSIONS 


This  report  describes  research  on  geodetic  algorithms, 
with  emphasis  on  the  prediction  of  gravity  anomalies  and  de¬ 
flections  of  the  vertical.  Two  general  types  of  methods  are 
considered : 

•  Classical  integral  formulas 

•  Modern  statistical  estimation  (least- 
squares  collocation). 

The  mathematical  details  of  the  major  algorithms  in  each  cate¬ 
gory  are  given  in  Appendix  A. 

The  present  work  builds  upon  and  complements  work 
done  by  the  IAG  Special  Study  Group  3.90  on  "Evaluation  of 
Local  Gravity  Field  Determination  Methods."  TASC  was  invited 
to  participate  in  this  Study  Group  because  of  our  development 
of  the  GEOFAST  algorithm.  Recent  work  by  the  IAG  Study  Group 
(Refs.  20  and  21),  which  is  based  mainly  on  analysis  of  gravity 
anomaly,  deflection  of  the  vertical,  and  terrain  height  data 
from  the  White  Sands,  NM  area,  is  summarized  in  Chapter  3. 
Chapter  4  of  this  report  describes  TASC’s  analyses  of  the  topo¬ 
graphic  contribution  to  deflections  of  the  vertical  in  the 
White  Sands  area.  This  report  also  includes  a  new  derivation 
of  the  convolutional  approximation  to  terrain  corrections  (Chap 
ter  5  and  Appendix  B).  Our  tests  of  the  GEOFAST  algorithm 
using  synthetic  marine  and  local  gravity  data  are  described  in 
Chapter  6. 
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This  final  chapter  is  divided  into  three  sections. 
Section  7.1  is  a  summary  of  our  key  contributions.  Section  7.2 
reiterates  the  major  conclusions,  and  Section  7.3  contains 
recommendations  and  suggestions  for  further  research. 


7.1  SUMMARY  OF  KEY  CONTRIBUTIONS 


The  key  contributions  of  TASC's  research  on  geodetic 
algorithms  are: 


•  Spectrum  and  cross-spectrum  analyses  of 
terrain  and  gravity  anomaly  profiles 
from  the  White  Sands,  NM  test  area  (Chap¬ 
ter  4,  Appendix  C) 

•  Analysis  of  correlations  between  topogra¬ 
phy  and  observed  deflections  of  the  verti¬ 
cal  in  the  White  Sands  area  (Chapter  4, 
Appendix  C) 

•  A  new  derivation  for  the  convolution 
approach  to  terrain  corrections  showing 
that  this  approach  is  actually  more  accu¬ 
rate  than  had  been  claimed  in  the  recent 
literature  (Chapter  5,  Appendix  B) 

•  Detailed  evaluation  and  quantification 
of  the  accuracy  of  the  GEOFAST  algorithm, 
based  on  synthetic  marine  and  local  grav¬ 
ity  data  (Chapter  6) 

•  A  summary  in  consistent  notation  of  the 
major  geodetic  algorithms  (Appendix  A) 

•  A  compact  summary  of  recent  results  by  the 
IAG  Special  Study  Group  3.90  (Chapter.3). 


1 

I 

I 


These  contributions  fall  into  the  two  general  areas  of:  anal¬ 
ysis  of  free-air  gravity,  deflections  of  the  vertical,  and 
topographic  effects  in  the  White  Sands  area;  and  detailed  com¬ 
parison  and  evaluation  of  geodetic  algorithms  including  the 
GEOFAST  algorithm. 
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7 . 2  CONCLUSIONS 


The  spectral  analysis  of  free-air  gravity  and  topo¬ 
graphic  effect  profiles  reflects  the  strongly  anisotropic  char¬ 
acteristics  of  the  White  Sands  area.  In  particular,  the  data 
profiles  were  divided  into  four  categories,  bas^d  on  geographic 
orientation  of  the  profile  segment  (north-south  or  east-west) 
and  type  of  terrain  (mountain  or  plateau).  The  main  conclu¬ 
sions  from  power  spectrum  and  coherence  analysis  are: 

•  Plateau  profile  segments  show  consistently 
lower  power  for  topographic  effect  than 
for  free-air  anomalies 

•  Mountain  segments  show  higher  power  when 
the  profile  orientation  is  north-south 

•  East-west  segments  show  local  peaks  in 
both  coherence  and  power  that  are  absent 
in  the  north-south  segments 

•  Short  wavelength  coherence  between  topo¬ 
graphic  effect  and  free-air  anomaly  is 
high  for  east-west  segments  and  low  for 
north-south  profile  segments. 

Deflections  of  the  vertical  in  the  White  Sands  are 
predicted  from  topography  alone  show  a  high  correlation  (80- 
90  percent)  with  observed  (astrogeodetic )  deflections,  as  illus¬ 
trated  by  means  of  scatter  diagrams.  The  predictions  account 
for  about  50  percent  of  the  rms  and  80  percent  of  the  mean  in 
the  observed  data.  Results  are  not  improved  substantially  by 
assuming  isostatic  compensation  or  different  densities  for 
mountain  and  plateau  regions . 

Numerical  experiments  with  the  GEOFAST  algorithm  show 

that : 
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•  GEOFAST  determines  deflections  of  the 
vertical  to  within  0.2  arc  sec  rms  from 
synthetic  marine  gravity  data 

•  This  accuracy  degrades  to  0.3  arc  sec  in 
the  presence  of  5  mgal  (rms)  noise,  and 
to  1  arc  sec  if  a  global  geopotential 
reference  model  is  not  used 

•  The  predicted  dependence  of  computation 
time  on  number  of  data  points  and  on 
bandwidth  is  achieved  in  practice.  A 
3600-point  data  set  is  computed  in  4  min 
with  nig=0  and  in  30  min  with  oig=4. 

•  An  accuracy  of  0.5  arc  sec  rms  is  achieved 
with  the  synthetic  local  gravity  data 

set  (5  km  sample  spacing)  if  a  reference 
field  including  wavelengths  longer  than 
100  km  is  removed. 


Each  of  the  four  algorithms  applied  to  predict  deflec¬ 
tions  of  the  vertical  in  the  White  Sands  area,  as  described  in 
the  1985  Report  (Ref.  21),  is  capable  of  approximately  1  arc  sec 
(rms)  accuracy  in  both  components.  Despite  the  demonstrated 
generality  and  flexibility  of  least-squares  collocation,  algo¬ 
rithms  based  on  classical  integral  formulas  remain  viable, 
especially  in  combination  with  FFT  implementations. 


7.3  RECOMMENDATIONS  AND  SUGGESTIONS  FOR  FURTHER  WORK 


These  are  the  principal  findings  and  recommendations 
of  this  study: 


•  Algorithms  based  on  both  modern  estima¬ 
tion  (collocation)  and  classical  integral 
formulations  have  important  contemporary 
applications 

•  For  a  given  application,  the  following 
considerations  apply  in  selecting  an 
appropriate  algorithm: 
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-  Where  availability  of  error  esti¬ 
mates  is  an  important  consideration, 
collocation  methods  should  be  used 

-  When  input  data  are  of  mixed  types, 
collocation  should  be  used 

-  If  statistical  optimality  of  the 
estimates  is  an  important  consid¬ 
eration,  collocation  is  indicated 

-  When  a  large  number  of  estimates 
are  required,  the  FFT  implementa¬ 
tion  of  the  classical  integral  formu¬ 
las  and  the  GEOFAST  FFT  implementation 
of  collocation  are  expected  to  be 

the  most  efficient  choices 

-  For  a  small  number  of  estimates, 
local  collocation  and  non-FFT  imple¬ 
mentations  of  classical  integral 
formulas  are  most  efficient 

•  Since  the  achieved  levels  of  accuracy 
for  algorithms  in  the  IAG  Study  reflect 
primarily  data  quality  and  availability, 
the  ultimate  accuracy  capability  of  the 
algorithms  must  be  evaluated  using  syn¬ 
thetic  data  that  include  a  realistically 
broad  spectrum  of  wavelengths 

•  Algorithms  should  be  applied  to  residual 
data  with  known  terrain  effects  and  global 
reference  fields  removed 

•  Remove-and-restore  techniques  which  incor¬ 
porate  the  above  removal ,  plus  a  later 
restoration  of  the  terrain  and  global 
long  wavelength  effects,  are  essential 

to  realize  the  full  performance  of  both 
classical  and  modern  geodetic  algorithms. 


Since  available  computer  power  is  still  a  limitation 
on  the  accuracy  and  utility  of  geodetic  algorithms  in  practice, 
we  suggest  consideration  of  the  implementation  and  practical 
use  of  new  computer  architectures  such  as  vector,  parallel, 
and  systolic  processors. 
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APPENDIX  A 
GEODETIC  ALGORITHMS 


A. 1  INTEGRAL  EVALUATION  METHODS 


The  purpose  of  this  appendix  is  to  review  briefly  the 
classical  integral  evaluation  methods  for  determining  the  undu¬ 
lation  of  the  geoid  and  the  deflections  of  the  vertical  from 
gravity  data.  These  methods  include  the  original  formulations 
of  Stokes  and  Vening  Meinesz  as  well  as  more  recent  developments 
and  extensions.  They  are  treated  in  the  standard  reference 
literature  (Ref.  1,  for  example). 


The  classical  formulations  assume  the  availability  of 
input  quantities  (measured  gravity)  on,  or  close  to,  the  physi¬ 
cal  surface  of  the  earth.  Since  the  variables  appearing  in 
the  classical  solution  are  linearized  small  quantities  which 
often  correspond  to  the  difference  in  shape  and  gravity  field 
between  the  ellipsoid  and  the  geoid ,  it  is  frequently  neces¬ 
sary  to  reduce  measured  gravity  to  the  geoid  (free-air  reduction) 
for  comparison  with  normal  gravity  on  the  ellipsoid.  The  inte¬ 
gral  formula  of  Stokes  relates  the  separation  between  geoid  and 
ellipsoid  at  a  point  (the  undulation  N)  to  the  gravity  anomaly 
on  the  entire  geoid: 


where 


N  =  4  s 


(<J> )  da 


(A. 1-1) 


R  is  an  assumed  mean  radius  for  the  earth 
(usually  taken  as  the  radius  of  a  sphere 
having  the  same  volume  as . the  ellipsoid) 
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Y  is  an  assumed  mean  value  for  gravity 

a  represents  the  surface  of  a  sphere  approxi¬ 
mating  the  ellipsoid 

Ag  is  the  free-air  anomaly 

4>  is  the  angular  separation,  on  the  surface  of 
the  sphere,  between  the  point  for  which  the 
solution  is  being  computed  and  the  element 
of  integration,  do 

S(i|»)  is  the  Stokes  function 


The  Stokes  function,  S(tJ>),  can  be  expressed  in  the  form 


S(«|»)  = 


2n+l 


n-1  n 


P  (cos  t|<)  (A. 1-2) 


n=2 


or  in  the  equivalent  finite  form 


S ( )  =  1  +  cosec  (^)  -6  sin  <!>  -5  COS  tjj 


•3  cos  41  £n[sin  (|-)  +  sin^  (|-)  ] 


(A. 1-3) 


Similarly,  the  Vening  Meinesz  formulas  relate  the  north-south 
(4)  and  east-west  (n)  components  of  the  deflection  of  the  ver¬ 
tical  at  a  point  on  the  geoid  to  worldwide  gravity  anomalies: 

i  =  ff Ag  (at)  C0S  “  da  (A. 1-4) 

a 

n  =  4ry  ffAg  (a|)  sin  “  da 

a 
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where 


j—  ,  the  derivative  of  the  Stokes  function,  is 
the  Vening  Meinesz  function 

a  is  the  azimuth  of  the  great-circle  arc  (of 
length  tjj )  joining  the  solution  point  and 
the  element  of  integration,  do 


An  explicit  expression  for  the  Vening  Meinesz  function, 

dS  .  .  .  „  .  ,  r 

-jjjj-  >  is  given  in  Eq .  A. 1-6: 


dS 

di|» 


cos(l) 

- ~ — : —  +  8  sin  -  6  cos  (£-) 

2  sin2(|)  2 


3[l-sin(|) ] 
sin  4» 


+  3  sin  4i  £n( sin(|) 


+  sin2(|)] 


(A. 1-6) 


Improvements  in  the  accuracy  of  the  basic  Stokes  and 
Vening  Meinesz  formulas  have  been  achieved  by  modifications 
that  take  into  account  the  effects  of: 

•  Gravitational  attraction  of  the  mass  of 
the  atmosphere  (Ref.  47) 

•  Local  topography  (Ref.  34) 

•  Integration  over  an  ellipsoid  instead  of 
the  assumed  sphere  (Ref.  48). 

Despite  these  improvements,  the  strongest  limitation  in  the 
use  of  these  formulas  is  that  they  require,  in  principle,  knowl¬ 
edge  of  gravity  anomalies  all  over  the  earth  in  order  to  compute 
undulation  or  deflection  at  a  particular  location.  Another 
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limitation  is  that  the  input  gravity  data  must  be  reduced  to 
the  geoid,  with  the  results  also  applying  to  this  theoretical 
surface  rather  than  to  the  actual  surface  of  the  earth.  Thus, 
for  example,  a  further  correction  for  the  curvature  of  the 
plumb  line  is  required  to  relate  the  results  of  the  Vening 
Meinesz  computation  to  deflections  actually  measured  at  the 
surface  of  the  earth. 


Approaches  to  addressing  these  limitations  will  be 
reviewed  below.  These  include  the  use  of  satellite-derived 
spherical  harmonic  data,  in  the  form  of  geopotential  spherical 
expansions,  to  replace  or  supplement  actual  gravity  anomaly 
data  in  remote  regions,  and  the  so-called  modern  or  nonclas- 
sical  approach,  introduced  by  Molodenskii,  that  uses  data  re¬ 
ferred  to  the  physical  surface  of  the  earth. 

Modern  (nonclassical )  solution  -  Associated  usually 
with  the  name  of  the  Russian  geodesist  Molodenskii,  the  so- 
called  modern  solution  to  the  problem  of  determining  the  shape 
of  the  earth  (and  hence  the  deflection  of  the  vertical)  from 
gravity  measured  on  the  surface  departs  from  the  classical 
approach  in  that  the  use  of  the  geoid  is  dispensed  with;  the 
desired  output  quantities  (such  as  height  anomaly  or  deflection 
at  the  surface)  as  well  as  the  input  quantities  (measured 
gravity)  are  referred  to  the  physical  surface  of  the  earth. 
The  modern  solution  was  introduced  in  1945  and  has  undergone 
considerable  development  and  improvement  since  that  time.  An 
up-to-date  reference  source  for  the  modern  solution  is  Part  D 
of  Ref.  2,  while  earlier  presentations  are  found  in  Ref.  3  and 
Chapter  8  of  Ref.  1. 


THE  ANALYTIC  SCIENCES  CORPORATION 


Molodenskii ' s  original  formula  for  the  determination 
of  the  deflection  of  the  vertical  from  worldwide  gravity  data 
takes  the  form  (Ref.  1): 

{  =  zs7  Jjf  <^+Gi)  (at)  cos  “  d°  -  f4  tan  Pi  <A-i-7> 

a 

H  =  5^  jy  [g|]  sin  a  do  -  tan  (A. 1-8) 

CT 

where  the  symbols  have  the  same  meaning  as  in  the  correspond¬ 
ing  Eqs.  A. 1-4  and  A. 1-5,  except  as  indicated  in  the  discussion 
to  follow. 

The  free-air  anomaly  Ag  appearing  in  Eqs.  A. 1-7  and 
A. 1-8  differs  from  the  classical  free-air  anomaly  (as  used  in 
Eqs.  A.  1-4  and  A.  1-5)  in  that  it  is  an  anomaly  referred  to 
ground  level  rather  than  to  the  geoid: 

Ag  =  gp  -  Yq  (A. 1-9) 

where 

g  is  actual  gravity  measured  at  a  point  P 
p  on  the  surface  of  the  earth 

Yq  is  normal  gravity  at  a  point  Q,  on  the 
^  reference  surface  known  as  the  telluroid 
(Ref.  1),  corresponding  to  P. 

Furthermore,  the  anomaly  referred  to  ground  level  is  modified 
by  ,  which  is  a  terrain  (and  local  mass  anomaly)  correction 
term  to  account  for  the  departure  of  the  earth's  physical  sur¬ 
face  from  a  surface  of  equal  potential.  This  term  is  calcu¬ 
lated  by  approximating  an  integral  of  the  form: 
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O  o 


which  involves  the  terrain  relief  (h-hp)  between  the  point  for 
which  the  correction  is  being  determined,  P,  and  other  points  in 
its  vicinity;  as  well  as  the  gravity  anomaly  Ag.  In  Eq.  A. 1-10, 
the  chordal  distance  between  the  point  P  and  the  element  of 
integration  is  the  other  symbols  have  been  defined  above. 

In  theory,  the  integral  in  Eq.  A.  1-10  is  extended  over  the 
entire  sphere.  Practical  implementations,  however,  generally 
do  not  extend  beyond  50  km  from  the  computation  point. 

The  second  term  in  Eqs .  A. 1-7  and  A. 1-8  represents 
the  effect  of  the  slope,  or  tilt,  of  the  terrain  in  the  vicinity 
of  the  computation  point,  with  p^  representing  the  inclination 
of  a  north-south  terrain  profile  with  respect  to  horizontal, 
and  p 2  representing  the  corresponding  east-west  inclination. 

The  Pellinen  Version  -  A  subsequent  modification  of 
the  Molodenskii  solution  by  Pellinen  (Ref.  4)  offers  computa¬ 
tional  advantages  as  well  as  a  simpler  formulation.  This  is 
the  formulation  which  is  the  basis  for  the  present  state  of  the 
art  of  gravimetric  determination  of  deflections  of  the  vertical 
using  a  deterministic  (nonstatistical )  approach. 

Pellinen' s  formulas 

4  =  sib  If (Ag)’  (at)  cos  “  da  (a- i'ii) 

a 

n  =  5^7  If  <Ag>'  (df)  sin  a  dCT  (A. 1-12) 

a 
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have  the  same  form  as  the  formulas  of  Vening  Meinesz  (Eqs.  A. 1-4 
and  A. 1-5),  except  for  the  use  of  the  modified  gravity  anomaly, 
(Ag) ' ,  in  place  of  the  classical  gravity  anomaly.  This  modified 
anomaly  incorporates  all  effects  of  the  topography  and  is  com¬ 
puted  as  follows: 


(Ag) '  =  Ag  +  C 


(A. 1-13) 


where 


Ag  is  the  free-air  anomaly  referred  to  ground 
level  (Eq.  A. 1-9) 

C  is  the  topographic  correction  term. 

The  correction  term  C  is  computed  in  theory  by  evaluating 
the  integral 

r2  Cf  (h-hp)(Ag-Agp) 

c  =  h  JJ — l-3 — ~d°  (A-1-U> 

a 

This  resembles  the  Molodenskii  terrain  correction  term,  G-^ ,  as 
defined  in  Eq.  A. 1-10,  except  for  the  use  of  a  relative  gravity 
anomaly ,  Ag-Agp,  referred  to  the  anomaly  at  the  computation 
point,  P.  In  practice  (as  was  pointed  out  in  the  discussion 
of  Eq.  A. 1-10),  the  integral  over  the  entire  sphere  appearing 
in  Eq.  A. 1-14  is  replaced  by  a  summation  limited  to  points  in 
the  immediate  vicinity  of  P. 

Numerical  implementation  -  The  practical  numerical 
implementation  of  formulas  like  Eq .  A. 1-1  for  the  determination 
of  undulation,  or  Eqs.  A. 1-11  and  A. 1-12  for  the  components  of 
deflection,  is  based  on  the  use  of  a  finite  sum  to  replace  the 
integral.  The  principle  will  now  be  illustrated  by  a  discussion 
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of  the  method  for  computing  deflection  of  the  vertical.  In 
the  Pellinen  approach  discussed  above,  there  are  two  sets  of 
integrals  that  must  be  evaluated  numerically.  At  each  gravity 
station  a  topographic  correction  term,  C,  is  computed  by  ap¬ 
proximating  the  integral  of  Eq.  A. 1-14,  using  gravity  and  ter¬ 
rain  relief  data  in  the  vicinity  of  the  gravity  station.  This 
is  applied  as  an  adjustment  to  the  ground  level  free-air  anomaly. 
With  the  gravity  data  thus  adjusted,  the  integrals  in  Eqs .  A. 1-11 
and  A. 1-12  are  then  evaluated  by  a  second  numerical  procedure, 
in  which  the  integral  over  the  entire  sphere  is  partitioned 
into  five  zones  at  increasing  radial  distances  from  the  computa¬ 
tion  point.  The  representation  and  treatment  of  the  gravity 
data,  as  well  as  the  details  of  approximating  the  integral  by 
a  finite  sum,  vary  from  zone  to  zone,  as  summarized  below. 

The  evaluation  of  both  types  of  integrals  enumerated 
above  is  usually  done  through  the  use  of  a  circular  template 
system  of  the  type  introduced  by  Rice  (Ref.  5).  In  its  origi¬ 
nal  form,  this  system  consisted  of  57  concentric  rings  sur¬ 
rounding  the  computation  poirrt,  the  innermost  ring  beginning 
at  a  distance  of  100  m  from  the  center,  and  the  outermost  ring 
ending  at  a  distance  of  1094. k  km  from  the  center  (see  Table  I 
of  Ref.  5).  Note  that  because  of  the  singularity  of  the  Vening 
Meinesz  function  when  ,  the  angular  separation,  approaches 
zero,  a  special  treatment  was  used  to  incorporate  the  effect 
of  gravity  stations  within  100  m  (see  the  discussion  of  the 
innermost  zone  below);  and,  in  his  original  treatment,  Rice 
did  not  account  for  the  effects  of  anomalies  beyond  1094.3  km, 
accepting  the  contribution  of  the  distant  zone  as  an  error 
source.  The  system  of  rings  is  further  divided  into  zones  by 
a  set  of  radial  lines  projecting  from  the  center  at  an  angular 
separation  of  10  deg.  There  are,  then,  a  total  of  2052  zones 
(excluding  the  innermost  zone)  in  the  Rice  template  system. 
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The  key  point  in  this  system  is  the  choice  of  the 
radii  of  the  circles  bounding  the  various  rings.  These  radii 
are  chosen  (based  on  the  variation  of  the  Vening  Meinesz  func¬ 
tion  with  distance  from  the  center)  so  that  each  zone  contributes 
the  same  amount  of  deflection  at  the  center  for  a  given  mean 
value  of  gravity  anomaly  over  the  zone.  In  particular,  each 
Rice  zone  contributes  0.001  arc  sec  of  radial  deflection  if 
the  mean  gravity  anomaly  for  the  zone  is  1.0  mgal . 

With  the  use  of  this  template  system,  the  numerical 
evaluation  of  an  integral  reduces  to  two  simple  steps: 

•  The  determination  of  a  mean  (or  other 
representative)  value  of  the  gravity 
anomaly  wiciiin  each  zone 

•  The  summation  of  the  contributions  from 
all  of  the  zones. 

For  the  evaluation  of  the  topographic  correction  term 
(Eq.  A. 1-14),  a  template  system  is  devised  in  Ref.  6  that  con¬ 
sists  of  90  zones,  having  the  characteristics  summarized  in 
Table  A. 1-1  (Ref.  6). 

Evaluation  of  the  topographic  correction  term,  using 
this  template  system,  is  a  two-stage  process: 

•  Determine  suitable  mean  values,  within 
each  zone ,  of 

gravity  anomaly  difference 
terrain  height  difference 

•  Compute  a  weighted  sum  of  the  individual 
zone  contributions. 

Reference  6  recommends  a  uniform  application  of  the  topographic 
correction  computation  to  each  gravity  station  in  the  point 
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TABLE  A. 1-1 


TEMPLATE  SYSTEM  FOR  TOPOGRAPHIC  CORRECTION 


RING  NUMBER 
(Note  1) 

CONTRIBUTION  TO 
TOPOGRAPHIC 
CORRECTION  (mgal) 
(Note  2) 

INNER 

RADIUS 

(km) 

OUTER 

RADIUS 

(km) 

NUMBER 

OF 

ZONES 

2 

7.5xl0-5 

0.100 

0.137 

6 

3 

7.5xl0‘5 

0.137 

0.217 

6 

4 

7.5xl0-5 

0.217 

0.526 

6 

5 

1.5xl0'5 

0.526 

0.734 

9 

6 

1.5xl0'5 

0.734 

1.218 

9 

7 

1.5xl0*5 

1.218 

3.556 

9 

8 

4.0xl0"6 

3.556 

7.287 

9 

9 

l.OxlO"6 

7.287 

9.878 

12 

10 

1 . Ox  10*  6 

9.878 

15.330 

12 

11 

l.OxlO'6 

15.330 

34.210 

12 

Note  1:  The  effects  of  Ring  1  (points  within  100  m) 
require  special  treatment  because  of  the 
divergence  of  the  integral 


Note  2:  This  is  the  contribution  to  the  correction 
term  from  a  zcne  for  which  the  mean  height 
difference  is  100  m  and  the  mean  gravity 
anomaly  difference  is  1  mgal . 


gravity  anomaly  data  base.  In  Ref.  7,  on  the  other  hand,  the 
correction  is  applied  much  more  selectively  (with  due  regard 
to  computational  resources)  at  points  where,  based  on  the  nature 
of  the  surrounding  terrain  and  the  amount  of  short-wavelength 
energy  in  the  gravity  anomaly  field,  the  topographic  correction 
is  expected,  a  priori,  to  be  significant. 
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The  first  step  in  the  evaluation  of  the  integral  formu¬ 
las  of  Eqs .  A. 1-11  and  A. 1-12  is  the  partitioning  of  the  domain 
of  integration  into  five  zones: 


•  Innermost  -  within  100  m  of  the  computa¬ 
tion  point 

•  Inner  -  a  3-deg  square  centered  on  the 
one-deg  square  containing  the  computation 
point 

•  Near  -  the  15-deg  square  centered  on  the 
5-deg  square  containing  the  computa¬ 
tion  point 

•  Middle  -  the  45-deg  square  centered  on 
(but  excluding)  the  near  zone 

•  Outer  -  the  rest  of  the  earth's  surface, 
beyond  the  middle  zone. 


Within  each  of  the  zones,  the  contribution  to  the 
integral  is  evaluated  as  outlined  below. 

Innermost  Zone  -  A  special  treatment  is  required  for 
the  innermost  zone,  which  must  be  excluded  from  the  domain  of 
integration  in  Eqs.  A. 1-11  and  A. 1-12  to  avoid  divergence  of 
the  integral,  since  the  Vening  Meinesz  function  behaves 
asymptotically  as 


dS 

di!> 


(A. 1-15) 


as  <f«  approaches  zero.  Within  this  innermost  zone  there  must 
be  sufficient  gravity  data  to  permit  an  estimate  of  the  north- 
south  and  east-west  components  of  the  horizontal  gradient  of 
the  gravity  anomaly.  The  traditional  approach  (Ref.  5)  is  to 
use  a  square  mesh  of  8  points  arranged  symmetrically  within 
the  innermost  zone  for  the  gradient  estimation.  Quite  clearly, 
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the  use  of  gradiometer  data,  if  available,  would  be  particularly 
appropriate  in  this  application.  The  innermost  zone  contribu¬ 
tions  are  then: 


^ IM  ~  "  2 y  (§rad  Ag)NS  rIM 


•  (A. 1-16) 


n 


A 


1 

2y 


(grad  Ag)EW  rIM 


(A. 1-17) 


where 


Y  is  normal  gravity  at  the  computation 

point 

rTM  is  the  radius  of  the  innermost  zone 

1  (typically,  r^  is  100  m) 

(grad  Ag)NS  is  the  north-south  component  of  the 
horizontal  gradient  of  the  gravity 
anomaly 


(grad  Ag ) Ey  is  the  east-west  component  of  the 
horizontal  gradient  of  the  gravity 
anomaly 


Inner  Zone  -  Within  the  inner  zone,  a  circular  template 
system  based  on  the  concept  of  Rice  rings  is  employed  to  organize 
the  computations.  The  smallest  of  the  21  rings  begins  at  a  dis¬ 
tance  of  100  m  from  the  computation  point,  while  the  outermost 
ring  ends  at  a  distance  of  128.6  km.  Since  the  rings  are  sub¬ 
divided  uniformly  by  radial  lines  separated  by  10  deg,  the  total 
number  of  compartments  is  756.  The  radii  defining  the  rings 
are  chosen  so  that  any  compartment  will  contribute  0.002  arc  sec 
of  radial  deflection  if  the  modified  gravity  anomaly  within 
that  compartment  has  a  uniform  value  of  1  mgal . 
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The  computational  procedure  for  the  inner  zone  consists 
of  these  stages: 


•  Sorting  the  point  gravity  anomalies  into 
the  appropriate  compartments 

•  Using  interpolation  procedures  to  compute 
a  best  estimate  of  the  gravity  anomaly 

at  the  center  of  each  compartment 

•  Summing  the  contributions  from  ail  of 
the  compartments. 


Beyond  the  inner  zone,  the  gravity  data  are  used  in 
the  form  of  mean  gravity  anomalies  defined  for  square  compart¬ 
ments.  The  contribution  from  the  area  of  overlap  between  the 
inner-zone  ring  template  and  the  block  pattern  of  the  near 
zone,  known  as  the  transition  area,  is  evaluated  by  an  elabo¬ 
rate  special  procedure  (Refs.  6  and  7)  to  ensure  that  no  area 
is  omitted  and  no  contribution  is  duplicated. 


Near  Zone  -  Within  the  near  zone,  the  gravity  field 
is  represented  by  0.5-deg  mean  gravity  anomalies;  the  integrals 
(Eqs.  A. 1-11  and  A. 1-12)  reduce  to  simple  sums  over  the  indi¬ 
vidual  blocks,  with  the  kernel  function 


or 


(A. 1-18) 


(A. 1-19) 


evaluated  at  the  center  of  the  block. 
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Middle  and  Outer  Zones  -  Mean  gravity  anomalies  de¬ 
fined  for  one-deg  squares  are  used  for  the  middle  zone;  five-deg 
mean  gravity  anomalies  are  used  in  the  outer  zone.  The  evalua¬ 
tion  of  the  contribution  of  each  individual  block  is  essentially 
as  described  for  the  near  zone. 

Summary  -  The  salient  features  of  the  evaluation  pro¬ 
cedure  described  above,  which  is  presented  in  detail  in  Refs.  6 
and  7,  are  summarized  in  Table  A. 1-2,  which  reviews  the  type 
of  data  and  the  computational  method  used  within  each  of  the 
zones . 


Use  of  spherical  harmonic  expansions  -  Although  the 
integral  formulas  for  undulation  and  deflection  nominally  re¬ 
quire  a  knowledge  of  point  gravity  anomaly  values  over  the 
entire  globe,  it  is  clear  that  the  contribution  from  zones 
distant  from  the  point  of  computation  comes  only  from  the  long- 
wavelength  content  of  the  gravity  field.  In  this  approach, 

TABLE  A. 1-2 

OVERVIEW  OF  THE  EVALUATION  OF  THE 
PELLINEN  INTEGRAL  FORMULA 


ZONE 

DATA  USED 

COMPUTATIONAL  METHOD 

- 

Innermost 

Point  gravity  anomalies 

Eqs.  A. 1-16  and  A. 1-17 

(note  possibility  of 

using  gradiometer  data) 

Inner 

Point  gravity  anomalies, 

Circular  template  system 

terrain  height'  data 

(Rice  rings) 

Near 

0.5-deg  mean  gravity 

Summation  over  indi- 

anomalies 

vidual  blocks 

Middle 

1-deg  mean  gravity 

Summation  over  indi- 

anomalies 

vidual  blocks 

Outer 

5-deg  mean  gravity 

Summation  over  indi- 

anomalies 

vidual  blocks 
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detailed  gravity  anomaly  data  are  used  only  within  an  annular 
cap  of  radius  4iq  centered  on  the  point  of  computation.  For 
the  computation  of  the  undulation,  for  example,  the  following 
approach  (Ref.  8)  is  representative: 

The  input  data  are  the  measured  surface  gravity  anoma 
lies  Ag  within  the  annular  cap  of  radius  (J/Q,  and  the  computed 
gravity  anomalies  Agg  outside  of  this  region.  The  computed 
anomaly  is  expressed  in  terms  of  the  spherical  harmonic  expan¬ 
sion  of  the  geopotential  in  the  form 

00  Q 

Ag  =  ^  V"1  (n-1)  (!)  Y''  (C*  cos  mA+S  sin  mA)P  (sin  $) 

6s  2  /  i  R  /  .  run  nm  run 

K  ~ 

n=2  m=o 

(A. 1-20) 

where 


GM  is  the  geocentric  gravitational  constant 

R  is  the  distance  of  the  computation  point 

from  the  earth's  center 


a  is  the  earth's  equatorial  radius 


$  and  A  are  geocentric  latitude  and  longitude 
of  the  computation  point 


P  is  the  associated  Legendre  function  of 

nm  .  .  , 

degree  n  and  order  m 


In  Eq.  A. 1-20,  S  and  C'  are  the  spherical  harmonic  coef- 

ficients  of  the  geopotential,  with  the  prime  signifying  that 

the  C  coefficients  have  been  adjusted  to  remove  the  effect 
run 

of  the  normal  gravity  field  associated  with  the  reference 
ellipsoid . 
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The  undulation  N  is  then  expressed  exactly  in  the  form 
N  =  Nx  +  N2  +  N3  (A. 1-21) 

where  is  the  imdulation  that  would  be  computed  from  spherical 
harmonic  data  only;  ^  and  represent  the  contribution  to 
the  undulation  of  Ag-Ag  in  the  annular  cap  and  the  rest  of 
the  earth,  respectively.  The  contribution  is  written  ex¬ 
plicitly  as 

«  n 

N1  =  i?  2]  (i)  H  (Cnmcos  mA+Snmsin  m^)Pnm(sin 

n=2  m=o 


(A. 1-22) 


where  y  is  normal  gravity  at  the  computation  point.  The  annular 
cap  contribution  has  the  same  form  as  Eq.  A. 1-1: 


N. 


R 

4ny 


(Ag-Ags)  S(t|>)  da 


(A. 1-23) 


except  that  the  adjusted  anomaly  (Ag-Agg)  is  used  and  the  inte¬ 
gration  extends  only  over  the  annular  cap.  Finally,  the  utility 
of  this  approach  depends  on  choosing  the  radius  of  the  annular 
cap  so  that  the  contribution  of  can  be  neglected.  Hence, 
the  formula  used  in  practice  is 

N  Nx  +  N2  (A. 1-24) 

The  error  contributed  by  neglecting  can  be  reduced 
by  modifying  the  Stokes  kernel  S(t|i),  as  discussed,  for  example, 
in  Refs.  9  and  10. 
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A.  2  LEAST-SQUARES  COLLOCATION 


Least-squares  collocation  is  a  modern  technique  for 
estimating  geodetic  quantities  from  observed  data.  It  has 
several  advantages  over  classical  integral  formulas: 


•  Collocation  allows  a  mix  of  different 
types  of  quantities  in  the  input  data 
and  different  types  may  likewise  be  esti¬ 
mated.  For  example,  the  classical  Vening 
Meinesz  integral  permits  computation  of 
deflections  of  the  vertical  from  measured 
gravity  anomalies.  However,  least-squares 
collocation  may  use  observations  of  anoma¬ 
lies  and  other  quantities,  such  as  astro- 
geodetic  deflections  or  gravity  gradients, 
to  estimate  simultaneously  deflections  and 
additional  quantities,  such  as  undulation 
of  the  geoid. 

•  No  special  processing  is  required  for 
irregularly-spaced  observation  or  esti¬ 
mation  points. 

•  The  effects  of  errors  in  the  input  data, 
including  lack  of  independence,  are  easily 
included.  Also,  with  very  modest  extra 
computation,  the  variances  and  covariances 
of  the  estimated  quantities  may  be  deter¬ 
mined  . 

•  In  its  most  general  form,  collocation 
includes  the  estimation  of  parameters 
governing  systematic  variations  in  the 
data,  such  as  terms  pertaining  to  the 
normal  potential  and  reference  system  or 
the  density  used  for  terrain  corrections. 


However,  least- squares  collocation  has  its  disadvan¬ 
tages  ,  too : 


•  It  requires  the  inversion  of  a  square 

matrix  with  dimension  equal  to  the  number 
of  data  values.  This  limits  the  speed 
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of  the  algorithm  and  the  number  of  data 
that  may  be  processed. 

•  A  model  must  be  provided,  which  is  usually 
assumed  to  be  stationary  and  isotropic, 
for  the  covariance  of  the  anomalous  poten¬ 
tial.  (From  such  a  model  other  covariances 
and  cross-covariances  are  derived.) 

This  appendix  is  organized  into  five  subsections,  as 
follows.  Appendix  A. 2.1  briefly  outlines  the  development  and 
background  of  least-squares  collocation.  Appendix  A. 2. 2  de¬ 
scribes  the  least-squares  interpolation  of  gravity  anomalies, 
which  is  a  special  case  of  collocation.  It  covers  much  of  the 
mathematical  development  essential  to  collocation.  The  covari¬ 
ance  and  cross-covariance  of  geodetic  quantities  of  different 
types  are  covered  in  Appendix  A. 2. 3,  in  terms  of  linear  operators 
applied  to  the  covariance  of  the  anomalous  potential.  In  Ap¬ 
pendix  A. 2. 4,  least-squares  interpolation  is  generalized  to 
least-squares  collocation  involving  quantities  of  mixed  type 
by  combining  the  results  of  the  preceding  two  sections.  A 
final  generalization,  to  include  the  simultaneous  estimation 
of  linear  parameters,  is  mentioned  briefly.  Appendix  A. 2. 5 
concludes  with  a  discussion  of  the  practical  application  of 
collocation . 

A. 2.1  Background 

Least- squares  collocation  was  introduced  in  geodesy 
in  the  late  1960s  by  Krarup  (Ref.  11)  and  Moritz  (Refs.  12 
and  13).  It  is  a  generalization  of  least- squares  interpolation 
to  include  an  arbitrary  mix  of  types  of  geodetic  quantities 
and  the  simultaneous  linear  estimation  of  parameters  governing 
systematic  behavior  of  the  field.  Heiskanen  and  Moritz  (Ref.  1, 
pp .  267-270)  discuss  least-squares  prediction  of  gravity  anoma¬ 
lies  (using  the  term  prediction  to  encompass  both  interpolation 
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and  extrapolation),  although  their  derivation  is  limited  to 
error-free  data.  The  more  general  derivation  given  in  the 
following  subsection  includes  the  effects  of  measurement  errors 
in  the  data.  Independently,  Gentry  and  Nash  (Ref.  14)  devel¬ 
oped  an  optimal  algorithm  for  the  estimation  of  deflections  of 
the  vertical  from  gravimetry.  Like  least-squares  interpolation, 
it  is  a  special  case  of  collocation.  Tscherning  (Ref.  15) 
published  a  general  computer  program  (GEOCOL)  for  least-squares 
collocation  and  has  updated  it  frequently  ever  since.  Tait 
(Ref.  16;  see  also  Appendix  A. 3  and  Chapter  6  of  this  report) 
describes  a  fast  Fourier-domain  implementation  of  collocation 
for  application  to  points  on  a  regular  grid. 

A  fairly  recent  and  very  thorough  discussion  of  the 
fundamentals  and  theory  of  collocation  is  given  by  Moritz 
(Ref.  2);  his  review  paper  (Ref.  17)  provides  a  more  compact 
summary.  Published  accounts  of  various  applications  of  col¬ 
location  are  abundant  (e.g.,  Refs.  18,  19,  20,  21). 

A  remark  about  the  term  "least-squares  collocation" 
itself  may  be  helpful.  "Least-squares"  emphasizes  the  statis¬ 
tical  nature  of  the  algorithm,  but  "collocation"  by  itself 
generally  refers  to  a  strictly  deterministic  technique  used  in 
the  numerical  solution  of  inhomogeneous  differential  equations 
(Ref.  22).  In  such  equations,  a  differential  operator  L  applied 
to  an  unknown  function  u  equals  a  forcing  function  f. 

Lu  =  f  (A. 2-1) 

The  strategy  of  collocation  is  to  obtain  a  numerical  approxima¬ 
tion  to  the  solution  u,  in  terras  of  a  linear  combination  of 
basis  functions,  that  minimizes  | | Lu - f | |  where  the  norm  ||  || 

depends  on  the  values  of  Lu  and  f  only  at  a  set  of  discrete 
points,  rather  than  over  a  continuous  interval.  Such  a  norm 
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has  obvious  advantages  in  terms  of  numerical  computation.  In 
spite  of  major  differences,  there  are  parallels  between  this 
collocation  technique  and  least-squares  collocation.  Geodetic 
quantities  are  known  only  at  discrete  points  and  "least-squares" 
implies  minimization  of  a  particular  norm  involving  values 
only  at  those  points.  Furthermore,  collocation  may  be  viewed 
as  primarily  an  attempt  to  estimate  the  anomalous  potential  , 
from  which  all  other  quantities  may  be  derived.  In  this  light, 
the  anomalous  potential  is  the  analogy  of  the  function  u;  the 
observations  are  analogous  to  the  values  of  f  on  the  discrete 
point  set,  which  should  equal  linear  functionals  of  the  anoma¬ 
lous  potential,  i.e.,  the  analogs  of  Lu  evaluated  at  the  dis¬ 
crete  points. 

A. 2. 2  Least-Squares  Interpolation 

Least-squares  interpolation  is  a  special  case  of  least- 
squares  collocation  in  which  only  a  single  type  of  quantity  is 
involved  and  parameter  adjustment  is  not  included.  Here  it  is 
assumed,  to  be  concrete,  that  gravity  anomaly  observations  are 
to  be  interpolated.  Only  the  covariance  function  for  gravity 
anomalies  is  needed,  but  otherwise,  the  mathematical  development 
proceeds  just  as  in  least- squares  collocation  (without  parameter 
adjustment ) . 

Least- squares  interpolation  and  collocation  employ 
the  statistical  concepts  of  expectation  and  covariance.  The 
gravity-anomaly  covariance  function,  a  function  of  two  points  P 
and  Q,  is  defined  by 


C(P,Q)  =  E  { Agp  dgQ}  (A. 2-2) 

where  E{  }  is  the  expected  value  operator.  Frequently,  it  is 
assumed  that  the  covariance  function  is  stationary  and  isotropic 


A-  20 


THE  ANALYTIC  SCIENCES  CORPORATION 


with  respect  to  the  horizontal  coordinates,  which  means  that 
the  horizontal  coordinates  of  P  and  Q  enter  into  C(P,Q)  only 
through  the  horizontal  distance  between  P  and  Q.  On  a  global 
scale  this  is  the  spherical  distance  4/  defined  in  Fig.  A. 2-1 
for  two  points  on  the  unit  sphere,  where  9  and  A  denote  colati¬ 
tude  and  longitude,  respectively.  According  to  a  basic  relation 
in  spherical  trigonometry, 

cos  4)  =  cos  9  cos  9'  +  sin  9  sin  9'  cos(A-A’)  (A. 2-3) 

If  the  covariance  function  is  stationary  and  isotropic,  and 
if  r  denotes  radial  distance  from  the  center  of  the  earth,  then 

C(P,Q)  =  average  of  Ag(r,9,A)  i\g(r '  ,d  '  )  over 

all  pairs  of  points  such  that  r  =  rp , 
r'  =  rQ,  and  <l>  =  4>pQ 

C(rp>  Tq  ,  4>pq)  (A.  2-4) 


A-21 


THE  ANALYTIC  SCIENCES  CORPORATION 


The  gravi ty- anomaly  covariance  function  is  also  defined 
on  a  local  scale.  In  the  stationary  isotropic  case  it  depends 
on  the  Euclidean  horizontal  distance 

s  =  /(x-x ' +  (y-y ' ) ^  (A. 2-5 ) 

where  x  and  y  denote  horizontal  Cartesian  coordinates.  Letting 
z  denote  the  elevation 

C(P,Q)  =  C(zp ,  Zq,  s  Pq )  (A. 2-6) 

Now  consider  observations  of  the  gravity  anomaly  at 
the  points  Qn .  The  spatial  distribution  of  the 

observation  points  is  unimportant  and  may  be  quite  irregular. 
Each  observation  £  will  consist  of  the  true  gravity  anomaly  Ag 
plus  a  measurement  error  n.  At  the  point  , 

SLi  =  Ag.j^  +  n^  (A. 2-7) 

Given  the  covariance  function  of  the  gravity  anomaly  and  a 
covariance  matrix  for  the  measurement  errors,  the  remainder  of 
this  section  demonstrates  how  to: 

•  Find  the  optimal  estimate  of  Ag  at  the 
point  P 

•  Determine  the  variance  of  the  estimate 

•  Determine  the  covariance  of  estimates  of 
Ag  at  points  P^,P2».-.»  Pm- 

The  estimate  of  Ag  at  P,  Agp,  is  assumed  to  be  d  linear 
combination  of  the  observations 

n 

Agp  =  ^  ai  £.  (A. 2-8) 

i  =  l 
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The  error  in  this  estimate  is  defined  as 


eD  -  ASV  ~ 


(A. 2-9) 


where  Agp  is  the  true  value  of  the  gravity  anomaly  at  P.  The 
optimal  estimate  Agp  is  obtained  by  choosing  the  coefficients 
ai,a2,...>  an  to  minimize  the  expectation  of  the  squared  error 


n  n 


E{e2}  =  E{ Agp  2Agp  Ev.'EEvav 

i=l  J=1 


i  =  l 


n  n  n 

CPP  "  2  2  3iCPi  +  E  E  aiaj 


i=l 


i=l  j=l 


-a  .C.  . 
“  '  iJ 


where 


(A. 2-10) 


Cpp  =  E{ Agp} 


(A. 2-11) 


Cpi  =  E{Agp£i}  =  E{AgpAgi)  +  E{Agpni} 


(A. 2-12) 


Cij  =  E^i£j*  =  E{Ag.Agj}  +  E{ Ag^n j } 


+  E{Agjni}  +  E{ninj} 


(A. 2-13) 


Defining  the  measurement  error  covariance  as 


E{ninj } 


(A. 2-14) 
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(which  is  a  diagonal  matrix  if  the  measurement  errors  are  un¬ 
correlated)  and  assuming  that  the  true  values  and  measurement 
errors  are  independent, 

Cpp  =  C(P,P)  (A. 2-15) 

Cpi  =  C(P,Q.)  (A. 2-16) 

eij  =  cu  +  Dij  <A-2-17) 

where 

Cij  =  C(Qi'Qj)  •  (A. 2-18) 

2 

Now,  for  i  =  1,2,...,  n,  the  derivative  of  E{ep}  with  respect 
to  a^  is  set  to  zero,  yielding 

n 

E  Vj  =  cpi  t=i -2 . "  (a-2-i9> 

j=i 

The  solution  of  the  linear  system  in  Eq .  A. 2-19  may  be  substi¬ 
tuted  back  into  Eq .  A. 2-8  to  yield  the  optimal  estimate  Agp 
and  into  Eq.  A. 2-10  to  determine  the  variance,  completing  two 
of  the  three  tasks  listed  earlier.  These  results  are  best 
expressed  explicitly  in  matrix  notation,  which  is  also  helpful 
in  performing  the  third  task  --  determining  the  covariance  of 
multiple  estimates  of  Ag. 
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The  vector 


and  matrices 

A  =  (ax  a2  . . .  an) 
CP SL  ~  (CP1  CP2  CPn)  =  E{AgpiT} 

and 

^SLSL  ~  j  ^  =  E^-  -  ^ 


(A. 2-20) 


(A. 2-21) 

(A. 2-22) 

(A. 2-23) 


are  needed.  Note  that  A  and  Cp£  are  1  by  n  row  vectors.  In 
the  case  of  simultaneous  estimation  at  m  different  locations, 
they  generalize  to  m  by  n  matrices.  Equations  A. 2-8,  A. 2-10, 
and  A. 2-19  are  equivalent  to 

Agp  =  M  (A. 2-24) 

E{ 6p}  =  Cpp  -  2ACp£T  +  AC££AT  (A. 2-25) 

and 

£itAT  =  cj2  (A. 2-26) 


respectively. 
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Substitution  of  the  formal  solution  of  Eq.  A. 2-26 


A  = 


(A. 2-21) 


into  Eqs.  A. 2-24  and  A. 2-25  yields  the  optimal  estimate 


(A. 2-28) 


and  its  variance 


E{ej] 


'PP 


(A. 2-29) 


When  gravity  anomalies  are  interpolated  at  multiple 
points  P^»P2,,,,»  Pm»  the  optimal  estimate  and  variance  of 
each  point  are  still  given  by  Eqs.  A. 2-28  and  A. 2-29,  but  in 
addition,  the  covariance  between  errors  in  prediction  at  dif¬ 
ferent  points  can  be  computed.  Defining  the  vector  of  esti¬ 
mated  gravity  anomalies  as 


Ug  at  Px 
[  4g  at  P2 


(A. 2-30) 


Y*  at  Pm  J 


and  applying  Eq.  A. 2-28  to  each  point  yields 


'sZ 


SL 


(A. 2-31) 
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where  is  an  m  by  n  matrix  defined  by 


CS2  =  {C(P.,Q.)}  =  E{ s  £  j  (A. 2-32) 


In  terms  of  £  and  the  vector  of  true  gravity  anomalies, 


Ag  at  P, 


Ag  at  P, 


Ag  at  PB 


(A. 2-33) 


the  vector  of  errors  at  each  point  P^,P2»...,  Pm  is 


e  =  s  -  s 


-  "  Cs£  ^££  - 


(A. 2-34) 


The  m  by  m  covariance  matrix  of  errors  in  the  estimated  anomalies 

T 

is  simply  the  expected  value  o fee.  This  may  be  expanded  as 


E{e  eT}  =  E{  s  sT}  -  E{s  £T}  C'J  cj 


Cs£  Cji  Efi  sT) 


+  C,t  E{£  iT}  C-J  Cj 


(A. 2-35) 
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Defining 


Css  =  E{s  sT}  =  fC(P.,P..)}  (A. 2-36) 


recalling  Eq .  A. 2-23  and  Eq .  A. 2-32,  and  noting  that 


EU  sT}  =  Cj 

yields  the  final  result 


(A. 2-37) 


E{  e  e1}  =  C 


ss 


-  C 


s2 


r'1  r  T 
s£ 


(A. 2-38) 


Note  that  the  diagonal  terms  of  Eq.  A. 2-38  are  the  variances 
of  the  estimates  of  each  point  and  are  equivalent  to  Eq.  A. 2-29, 
while  the  off-diagonal  terms  are  covariances  between  estimates 
at  different  points. 

A. 2. 3  Covariance  Functions 

In  generalizing  least- squares  interpolation  to  least- 
squares  collocation  it  is  necessary  to  determine  the  covariances 
and  cross-covariances  of  quantities  of  different  types  (e.g., 
the  cross-covariance  between  north  deflection  of  the  vertical 
4  and  gravity  anomaly  kg).  K(P,Q),  the  covariance  function  of 
the  anomalous  potential  T,  is  central  to  this  process;  its  use 
guarantees  a  consistent  set  of  covariances  and  cross-covariances . 
In  analogy  to  Eq.  A. 2-2, 

K(P,Q)  =  E{TpTp}  (A. 2-39) 
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The  purpose  of  this  section  is  to  explain  the  derivation  of 
additional  covariances  from  K(P,Q)  and  to  present  some  examples 
of  stationary  isotropic  models  for  K(P,Q). 

Observable  quantities  such  as  gravity  anomaly,  geoid 
undulation,  deflection  of  the  vertical,  and  components  of  the 
gravity  gradient  depend  on  the  anomalous  potential  through 
linear  operators.  For  example,  the  relation  (Ref.  1) 

Ag  =  -  §|  -  §  T  (A. 2-40) 

may  be  expressed  as 

Ag  =  LAg  T 

where  the  notation  LAg  denotes  the  linear  operator: 

LAg  =  -  ~  -  |  (A. 2-41) 

Evaluation  at  a  specific  point  P  may  be  denoted  by  adding  a 
subscript  to  the  operator,  e.g., 

Agp  =  LpS  T  =  ('  I?  '  ?)  T(r,0,\)  (A. 2-42) 

r_rP 
8  “6p 

\=\p 

Geodetic  quantities  and  their  corresponding  operators  are  sum¬ 
marized  in  Table  A. 2-1.  (y  denotes  normal  gravity.) 
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TABLE  A. 2-1 

GEODETIC  QUANTITIES  AND  LINEAR  OPERATORS 


QUANTITY 

SYMBOL 

OPERATOR 

Anomalous  Potential 

T 

1  (the  identity) 

Gravity  Anomaly 

Ag 

3 _ 2 

3  r  r 

Geoid  Undulation 

N 

1/Yo 

|  north 

i 

1  3 

Y  R  36 

1  o 

Deflection  of  / 

the  Vertical  \ 

1  A 

1  east 

n 

X  u 

V 

Y  r  sin  0  3\ 
o 

Gradient  Components 

-- 

2  2 

3  /3r  ,  etc 

The  operators  in  Table  A. 2-1,  which  yield  geodetic 
quantities  when  applied  to  the  anomalous  potential,  may  be 
applied  in  pairs  to  the  covariance  function  of  the  anomalous 
potential  to  yield  covariances  and  cross-covariances  as  needed 
This  fundamental  observation  may  be  demonstrated  with  an  ex¬ 
ample.  Consider  the  cross-covariance  between  the  gravity 
anomaly  Ag  at  point  P  and  the  deflection  of  the  vertical  %  at 
point  Q, 

E{Agp4Q}  =  E{ ( Lp®  T)(l|  T>}  (A. 2-43) 

The  expectation  operator  is  linear  (it  is  a  kind  of  average), 
so  the  linear  operators  Lp®  and  Lq  may  be  moved  to  the  left. 

Then,  from  the  definition  of  K(P,Q)  (Eq.  A. 2-39), 
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E{AgpCQl  =  LpS  K 


(A. 2-44) 


where  Lpg  operates  on  K  as  a  function  of  the  coordinates  of  P 

and  Lq  operates  on  K  as  a  function  of  the  coordinates  of  Q. 
For  example,  the  above  cross-covariance  is  written  out  as 


LpSlqK=  <9VaF>  K(r,e,X,r',0,,V) 


r=rp  >  9=8p ,  A.=A-p 
r ' =rQ , 0 ' =0Q , A ' =\Q 


(A. 2-45) 


Another  example  is  the  relation  between  the  gravity  anomaly 
covariance  function  (Eq.  A. 2-2)  and  the  covariance  of  the 
anomalous  potential , 


C(P,Q)  =  Lpg  Lq8  K  (A. 2-46) 

In  practice,  K(P,Q)  is  usually  approximated  by  a  para¬ 
metric  model.  The  definition  of  K(P,Q),  Eq .  A. 2-39),  imposes 
some  constraints  on  such  models. 

K(P,Q)  =  E{ TpTg]  =  E{ TqTp}  =  K(Q,P)  (A. 2-47) 

so  the  covariance  function  must  be  symmetric  in  its  dependence 
on  the  coordinates  of  P  and  Q.  Also,  the  harmonicity  of  T  (it 
satisfies  Laplace's  equation  outside  the  earth)  implies  that  K 
must  be  harmonic  in  both  the  coordinates  of  P  and  the  coordinates 
of  Q.  Parametric  models  for  K  are  almost  always  taken  to  be 
stationary  and  isotropic,  and  in  addition,  they  must  be  positive 
semi-definite  to  ensure  positive  semi -de f ini te  covariance 
matrices . 
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All  global  covariance  functions  of  the  form 


Lr  n  ^ 

k  (— £-)  Pn  (cos  4> ) 

n  rPrQ  n 


(A. 2-48) 


where 


COS 


is  the  radius  of  the  earth, 

is  the  Legendre  polynomial  of  degree  n, 

is  given  by  Eq.  A. 2-3)  with  0p,  Ap,  0Q, 
Aa  substituted  for  0,A,8',A',  1 


are  stationary  and  isotropic,  because  the  horizontal  coordinates 
enter  in  only  through  the  distance  4> .  In  addition  they  are 
harmonic,  because,  for  all  non-negative  n, 


Pn  (cos  *> 


satisfies  Laplace's  equation  in  the  coordinates  of  both  P  and  Q, 
and  the  symmetry  between  P  and  Q  is  evident.  Such  covariance 
functions  are  positive  semi-definite  if  and  only  if 


k  >  0 
n  - 


(A. 2-49) 


Two  examples  of  covariance  functions  in  this  form  are 
the  Tscherning  and  Rapp  model  (Ref.  23)  in  which 


.  _  AR2Sn+1 

n  ~  (n-1 ) (n-2) (n+24)  ’ 


n  >  2 


(A. 2-50) 
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and  the  attenuated  white  noise  model  of  Heller  and  Jordan 
(Ref.  24),  where 


__  (2n+l)D2(2R-D)oT2  ,R.D.n+l 

n  (R-D)2(2R2-2DR+D2>  '  R  ) 


The  Tscherning  and  Rapp  model  has  two  parameters,  A 
and  S.  A  sets  the  overall  level  of  the  covariance  and  usually 
S  is  slightly  less  than  unity.  Also,  is  set  to  zero  for  n 
less  than  2,  and  for  very  large  values  of  n  as  well. 

The  attenuated  white  noise  model  also  has  two  param¬ 
eters.  D  is  the  depth  of  a  sphere  on  which  the  anomalous 

2 

potential  is  a  white  noise  process  and  aT  is  the  variance  on 
that  sphere.  In  addition  to  this  clear  physical  interpretation, 
the  covariance  function  for  the  attenuated  white  noise  model 
has  the  advantage  of  a  closed- form  expression,  which  is  more 
efficient  for  computational  purposes  than  the  summation  in 
Eq.  A. 2-48. 


Given  the  covariance  function  K(P,Q)  in  the  form  of 
Eq .  A. 2-48,  the  gravity  anomaly  covariance  may  be  expressed  in 
a  related  series.  It  is 


C  ( P  , 0)  = 


Z 


(n-l )  k 


R 


_Rjn+2  P 

rPrQ  ) 


(cos  4i) 


(A. 2-52) 


which  follows  from  Eq.  A. 2-46). 

For  local  covariance  functions,  in  the  flat-earth 
approximation,  there  are  analogs  to  the  above  expressions. 
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K(P 


as 

,Q)  =  / 


-W(2p+Z0) 

k(ui)e  ^  Jq(u)s)  u)du) 


(A. 2-53) 


where  J  is  the  zero-order  Bessel  function  of  the  first  kind  and 
o 


j. 

((Xp-xQ)2  +  (yp-yq)2] 


(A. 2-54) 


is  clearly  stationary,  isotropic,  and  symmetric.  It  is  also 
harmonic  in  the  coordinates  of  both  P  and  Q  because 


-uj(zp+zn) 
e  ^  JQ(u>s) 


satisfies  Laplace's  equation,  and  it  is  positive  semi-definite  if 
and  only  if 


k(u> )  >  0  (A. 2-55) 

for  all  w .  The  corresponding  expression  for  the  gravity  anomaly 
covariance  function  is 


C(P,Q) 


2  -w(zp+z  ) 

uj  k  (  uj  )  e 


J  (ius)  tuduj 
o 


(A. 2-56) 


Two  examples  are  the  flat-earth  form  of  the  attenuated 
white  noise  model  (Ref.  24)  in  which 

k(w)  -  4D  a^,  e  (A. 2-57) 


A-  34 


THE  ANALYTIC  SCIENCES  CORPORATION 


and  Reilly's  covariance  model  (Refs.  25,  26),  in  which 


k(w) 


CQ(d4/2)  e‘w2d2/2 


(A. 2-58) 


A. 2. 4  Collocation 


Given  a  covariance  function  for  the  anomalous  poten¬ 
tial  K(P,Q)  and  the  linear  operators  in  Table  A. 2-1,  it  is 
straightforward  to  generalize  least-squares  collocation  to 
estimate  quantities  of  diffe  nt  types  from  data  of  different 
types.  In  fact,  Eqs .  A.?  ,j.  and  A.  2-38  hold  exactly  if  the 
definitions  of  the  obser  ation  vector  .£,  the  estimation  vector 
s,  and  the  covariance  matrices  are  modified.  This  is  most 
easily  demonstrated  with  a  concrete  example. 


Suppose  that,  rather  than  consisting  solely  of  gravity 
anomaly  measurements  as  in  Eq .  A. 2-7,  the  observation  vector  £ 
includes  deflections  of  the  vertical  also,  e.g., 


at 

9l 

+ 

n 

at 

^2 

+ 

n 

\  :  / 


(A. 2-59) 


and  that  both  deflection  of  the  vertical  and  geoid  undulation- 
are  to  be  estimated,  e.g., 


1  i  at  P1^ 
N  at 

:  I 


(A. 2-60) 
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A. 2-31  and 


(A. 2-62) 


(A. 2-63) 


THE  ANALYTIC  SCIENCES  CORPORATION 


where  s 


I  true  value  of  £  at  P^\ 
true  value  of  N  at  P2 

:  I 


(A. 2-64) 


T 

and  E{ (s-s)  (s-s)}  has  been  minimized. 

By  the  above  approach,  the  least-squares  interpolation 
of  a  set  of  data  (including  errors)  of  homogeneous  type  has 
been  extended  to  least-squares  collocation,  in  which  the  types 
of  data  and  estimated  quantities  may  be  heterogeneous.  In 
addition,  it  is  possible  to  incorporate  the  simultaneous  esti¬ 
mation  of  a  set  of  parameters  upon  which  the  data  are  linearly 
dependent.  Such  systematic  adjustments  might  pertain,  for 
example,  to  the  geodetic  reference  system,  the  density  chosen 
for  terrain  corrections,  or  arbitrary  long-wavelength  trends 
in  the  data.  In  practice,  however,  such  adjustments  are  often 
done  independently  of  the  collocation  procedure  (e.g.,  Refs.  20, 
21)  and  the  details  are  omitted  here.  They  are  available  in 
many  of  the  works  by  Moritz  (e.g.,  Refs.  2,  17). 

A. 2. 5  Practical  Considerations 

A  remark  on  the  computational  aspect  of  least-squares 
interpolation  and  collocation  is  necessary.  The  matrix  inverse 


appears  in  both  the  expression  for  the  optimal  value  of  Agp 
(Eq.  A. 2-28)  and  the  expression  for  its  covariance  (Eq.  A. 2-29), 
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and  also  in  the  more  general  expressions  in  Eqs .  A. 2-31  and 
A. 2-38.  However,  in  practice,  it  is  more  efficient  and  accu¬ 
rate  to  avoid  the  explicit  computation  of  a  matrix  inverse. 
Instead,  one  should  deal  directly  with  the  fundamental  expres¬ 
sions  in  Eqs.  A. 2-24,  -25,  and  -26. 

Because  C££  is  symmetric  and  positive  semi-definite 
it  has  a  Cholesky  decomposition  (e.g,,  Ref.  27)  such  that 

C££  =  UTU  (A. 2-65) 

where  U  is  an  upper  triangular  matrix.  Once  U  has  been  com¬ 
puted,  Eq.  A. 2-26  is  solved  for  the  1  by  n  coefficient  matrix  A 
in  two  steps.  First,  the  solution  y  to  the  system 

UTy  =  cj£  (A. 2-66) 

is  computed.  This  involves  only  straightforward  substitution 

T 

because  of  the  triangular  structure  of  U  .  Next,  A  is  computed 
from 


UAT  =  y  (A. 2-67) 

which  is  likewise  straightforward.  Then  Agp  is  computed  directly 

from  Eq.  A. 2-24  rather  than  from  Eq.  A. 2-28.  Similarly,  E{ep} 
is  computed  from 

E{ Sp}  =  Cpp  -  AcJ£  (A. 2-68) 

(obtained  by  using  Eq .  A. 2-26  to  eliminate  C££  from  Eq .  A. 2-25) 
rather  than  from  Eq.  A. 2-29.  This  procedure  is  readily  extended 
to  compute  estimates  of  multiple  quantities  and  their  covariances 
(in  either  interpolation  or  collocation),  denoted  formally  by 
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Eqs.  A. 2-31  and  -38.  Note  that  the  Cholesky  decomposition  need 
not  be  repeated  because  the  same  data  covariance  matrix  1S 

needed  for  each  of  the  estimates. 

A  second  remark  concerns  the  covariance  of  the  gravity 
field  and  the  handling  of  data.  As  mentioned  above,  stationary, 
isotropic  covariance  models  are  normally  used.  However,  the 
gravity  field  of  the  earth  is  far  from  stationary.  Its  character 
varies  with  location  on  a  global  scale  (e.g.,  between  ocean 
basin,  continental  shelf,  continent,  etc)  and  on  a  local  scale 
''e.g.,  between  a  mountain  and  adjacent  lowland).  Also,  the 
prominence  of  linear  features  such  as  ocean  trenches  and  rises, 
island  chains,  mountain  chains,  and  continental  rift  zones 
indicates  that  the  field  is  far  from  isotropic.  In  practice 
(Ref.  21),  this  difficulty  is  circumvented  not  by  using  non¬ 
stationary,  anisotropic  covariance  models  (which  would  require 
too  many  parameters  to  be  practical)  but  rather  by  estimating 
residual  quantities  from  residual  data.  The  residuals  are 
computed  by  removing  systematic  components,  such  as  those  com¬ 
puted  from  topography  and  from  long-wavelength  models  of  the 
gravity  field,  from  the  data.  Collocation  is  performed  on  the 
residual  field,  which  is  generally  more  stationary  and  isotropic 
than  the  field  itself.  After  collocation,  the  systematic  com¬ 
ponents  are  restored,  i.e.,  added  to  the  estimated  residual 
quantities . 


A. 3  THE  GEOFAST  ALGORITHM 

A. 3.1  Formulation  of  the  Problem 

The  GEOFAST  algorithm  (Ref.  16)  was  developed  to  pro¬ 
vide  an  efficient  computational  solution  to  the  minimum- 
variance  estimation  equations 
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In  geodetic  work,  this  formula  is  commonly  known  as  least-squares 

collocation.  It  is  equivalent  to  Eq.  A. 2-31  with  minor  nota- 

tional  changes:  x  replaces  s,  z  replaces  i,  and  the  bar  on 

Czz  is  dropped  for  simplicity.  Here,  z  is  a  data  vector  of 

dimension  N  and  C  is  the  sum  of  its  auto-covariance  matrix 

zz 

(of  dimension  N  by  N)  plus  a  measurement  noise  covariance  matrix. 
The  vector  of  estimates  x  is  of  dimension  M,  and  Cxz  is  the 
cross-covariance  matrix  between  x  and  z  of  dimension  M  by  N . 

The  data  covariance  is  generally  modeled  by  an  analytic  function 
with  parameters  fitted  to  global  or  local  data.  The  cross¬ 
covariance  matrix  C  is  then  analytically  derived  from  the 
auto-covariance  by  using  the  physical  relationship  between  the 
estimated  quantities  x  and  the  observed  quantities  z  (Ap¬ 
pendix  A. 2. 3).  The  method  is  implicit  since  this  relationship 
does  not  appear  explicitly  via  a  measurement  equation  of  the 
form  z  =  H  x.  Formally,  the  covariances  are  defined  by 

CZ2  =  Eli  zT)  ,  Cxz  =  El*  zT)  (A. 3-2) 

(equivalent  to  Eqs .  A. 2-23  and  A. 2-32)  where  E  is  ^he  expectation 
operator  and  both  vectors  are  assumed  to  have  zerj  means. 

Computer  solution  of  Eq.  A. 3-1  by  standard  techniques 

3 

involves  on  the  order  of  N  operations  (such  as  multiplications 
and  additions)  to  obtain  the  intermediate  vector 


and  on  the  order  of  MN  operations  to  yield  the  estimates 

x=  Cxz  u  (A. 3-4) 
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This  limits  the  feasible  application  of  these  methods  to  data 
sets  of  a  few  hundred  points  at  most,  whereas  thousands  of 
measurements  are  often  available  from  sensors  such  as  satel¬ 
lite  altimeters. 

With  some  assumptions  about  the  data  and  the  use  of 
frequency  domain  techniques ,  this  workload  can  be  reduced 
dramatically  to  the  order  of  N  log2  N.  These  assumptions  are 
as  follows.  First,  let  the  data  points  z(0,<|>)  be  given  on  a 
rectangular  grid  in  the  0  , 4>  plane  (Fig.  A. 3-1).  Second,  as¬ 
sume  that  the  covariance  between  two  data  points  is  a  function 
of  relative  position  only.  That  is,  the  covariance  function 
is  shift-invariant  (or  a  displacement  kernel) 

COV  (z(01,01),  z(02,$2)J  =  f2z(92“ei>  ^2"^1 }  (A- 3-5) 


ft-49247 


COV  (P,Q)  =  COV  (P',Q')  =  fZ2  (d2-dv  ) 
Figure  A. 3-1  Data  and  Covariance  Structure 
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Note  that  this  includes,  but  is  not  restricted  to,  the  case  of 
isotropic  models.  Third,  similar  conditions  are  assumed  for 
the  cross-covariance  and  estimate  points  x(8,$),  namely, 

COV  [x(e1 ,4>1 )  ,  z(e2,$2)]  =  fxz(02“0i»  *2'*1*  (A. 3-6) 

Fourth,  the  estimates  are  to  be  obtained  on  the  same  grid  as 
the  data. 

The  foregoing  assumptions  are  consistent  with  practical 
applications.  The  GEOFAST  algorithm  maps  data  given  on  a  limited 
region  of  a  spherical  (or  ellipsoidal)  surface  onto  a  plane 
with  very  little  distance  distortion.  The  data  are  then  gridded. 
Shift  invariance  then  corresponds  to  the  common  model  of  sta¬ 
tionary  statistics.  The  estimate  grid  may  actually  be  any 
translation  (that  is  any  offset  in  altitude,  latitude,  and 
longitude)  of  the  data  grid.  With  a  fast  algorithm,  it  is 
often  cheaper  to  compute  estimates  at  all  grid  locations  even 
if  only  a  subset  is  wanted. 

The  estimation  equation  (Eq.  A. 3-1)  takes  on  a  special 
structure  which  is  exploited  by  the  GEOFAST  algorithm.  The 
natural  form  for  the  data  is  an  n^  by  n2  matrix 

1  =  tzjlJ  »  zjk  =  z(0Q+jA0,  0o  +  kA4>)  (A. 3-7) 

where  j=0 ,1 , . . . ,n^-l  and  k=0 , 1 , . . . ,n2~l  corresponding  to  the 

data  grid  (Fig.  A. 3-1)  with  origin  (0q,0q;  and  mesh  size  (A0  , 

A0 ) .  The  total  number  of  grid  points  is  N  =  n^n2>  The  data 

matrix  Z  is  put  in  vector  form  z  by  listing  its  elements  in 

row-by-row  order.  With  this  convention  the  N  by  N  covariance 

matrix  C  assumes  a  block  structure  (see  Ref.  16) 
zz 
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C 

zz 


.  .  T 


nr1 


.  .  T 


nl"2 


(A. 3-8) 


-(n1-l) 


■-<nr2) 


in  which  there  are  n-^  by  blocks  and  each  block  is  of  dimen¬ 
sion  n2  by  n2 •  The  block  matrix  is  the  cross-covariance 
between  the  elements  in  two  rows  of  Z  which  are  k  rows  apart. 
This  depends  only  on  k  by  shift-invariance.  By  the  same  argu¬ 
ment  each  block  matrix  T^  is  of  the  form 


C0 


t 


-1 


t  O 

n2-2 


(A. 3-9) 


t-(n2-l)  t-(n2-2)  *  •  t0 


where  t  is  evaluated  using  Eq.  A. 3-5  or  A. 3-6. 

The  matrix  C zz  is  said  to  be  block  Toeplitz  with 
Toeplitz  blocks,  or  simply  block  Toeplitz.  Note  that  Czz  is 
symmetric  and  is  completely  determined  by  its  first  block  row 
(or  block  column),  whereas  the  blocks  (except  Tq )  are  not  sym¬ 
metric  and  are  determined  by  their  first  row  and  column.  Alto 
gether  the  definition  of  Czz  requires  fewer  than  2n^n2  =  2N  pa 

rameters.  The  structure  of  C  is  also  block  Toeplitz,  but 

xz 

symmetry  does  not  hold  in  general;  approximately  4N  parameters 
are  needed  to  define  it.  In  summary,  the  geodetic  estimation 
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problem  becomes  equivalent  to  the  inversion  and  multiplication 
of  block  Toeplitz  matrices. 


A. 3. 2  Outline  of  the  Algorithm 


The  GEOFAST  algorithm  achieves  its  efficiency  by  trans¬ 
forming  the  estimation  equations  into  the  frequency  domain 
where  an  accurate  approximation  can  be  made  to  reduce  the  work¬ 
load.  The  transformation  is  accomplished  with  the  Fast  Fourier 
Transform  (FFT),  which  requires  on  the  order  of  N  log2  N  opera¬ 
tions.  The  transformed  covariance  matrices  are  closely  approxi¬ 
mated  by  (block)  banded  matrices  resulting  in  a  solution  workload 
of  order  N.  A  simple  trade-off  between  approximation  accuracy 
and  computer  workload  is  controllable  by  choice  of  the  bandwidth. 


A  simplified  flow  chart  of  the  algorithm  is  shown  in 

Fig.  A. 3-2.  Input  to  the  computation  consists  of  the  n^  by  n2 

data  matrix  Z,  together  with  the  definition  of  the  covariance 

matrices  C  ,  C  .  As  seen  in  Eqs .  A. 3-8  and  A. 3-9,  each  co- 
xz  zz 

variance  can  be  specified  by  ^.nin2  parameters  which  are  natur¬ 
ally  arranged  as  a  2n^  by  2n2  correlation  matrix.  Specifically, 
let  the  Toeplitz  matrix  T  in  Eq.  A. 3-9  be  represented  by  the 
column  vector 


t 


t-(n2-l) ' 


(A. 3-10 ) 


of  dimension  2n2>  Then  each  block  matrix  T^  in  Eq .  A. 3-8  can 
be  represented  by  a  vector  t^,  and  the  matrix  C-zz  by  the  2n^ 
by  2n2  matrix 


T  =  0  t  ,  .  i 

zz  -  i — (n  -1)  i 

i  1  i 


1  t-  I  r  I  f  1 

•  • 1  1  i  i  Li  i 

i  — 1  i  -0  i  -1  t 

i  iii 


'  !  -n,-l 


(A. 3-11) 
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•  TRANSFORM  INPUT  I  ■  COMPUTE  FOURIER  ESTIMATES  1  ■  TRANSFORM  OUTPUT 


L 


FREQUENCY  DOMAIN 


J 


Figure  A. 3-2  Two-Dimensional  GEOFAST  Algorithm 


For  the  symmetric  case  (C  )  only  half  of  the  matrix  T  is 

4iA  ZZ 

required  with  t_k  =  t_k.  Total  computer  storage  for  input  is 
of  order  N  =  n^^. 

The  first  stage  of  the  algorithm  is  the  transformation 
of  the  inputs  to  the  frequency  domain.  The  FFT  is  simply  an 
efficient  implementation  of  the  complex  Discrete  Fourier  Trans¬ 
form  (DFT)  which  is  defined  in  two  dimensions  by 


z 


pq 


1 

Vn1n2 


zjk  exp 


-2ni 


(A. 3-12) 


Since  the  estimation  equations,  Eq .  A. 3-1,  are 
data  vector  z  it  is  convenient  to  express  Eq. 
equivalent  matrix  form, 


in  terms  of  the 
A. 3-12  in  the 


z  1  =  F  z 


(A. 3-13) 
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The  N  by  N  complex  two-dimensional  DFT  matrix  F  has  a  special 
structure  which  is  discussed  further  in  Ref.  16.  The  vector 
z'  in  Eq.  A.  3-13  is  the  complex  representation  of  the  data 
vector  z  in  a  finite  Fourier  series  of  complex  exponentials. 
Since  real  quantities  are  more  convenient  to  deal  with  in  a 
computer,  a  further  transformation  to  a  sine  and  cosine  series 
is  performed  by  combining  real  and  imaginary  parts.  This  cor¬ 
responds  to  the  matrix  operations  (redefining  z’ ) 

z'  =  H  F  z  (A. 3-14) 

where  H  is  a  very  sparse  complex  N  by  N  matrix.  The  product 
HF  is,  of  course,  a  real  matrix. 

Finally,  the  accuracy  of  the  banded  approximation, 
which  is  essential  for  the  success  of  the  algorithm,  is  criti¬ 
cally  enhanced  by  the  introduction  of  a  data  window.  In  matrix 
form  this  becomes  (redefining  z'  again) 

z'  -  H  F  W  z  (A. 3-15) 

where  W  is  a  real  N  by  N  diagonal  matrix.  That  is,  each  data 
value  z.,  is  multiplied  by  a  scalar  weight  w.,  before  applying 

J  K  J  K 

the  DFT.  This  windowing  compensates  for  the  finite  extent  of 
the  data  grid  in  a  manner  explained  in  Ref.  16.  A  representa¬ 
tive  window  function  in  two  dimensions  would  be  a  gaussian 
probability  density  which  is  peaked  at  the  center  and  tapered 
toward  the  edges  of  the  data  grid.  In  GEOFAST  an  optimal  one¬ 
dimensional  window  w°  derived  by  Kaiser  (Ref.  28)  is  used  to 
form  a  product  window  w ^ ^  =  (w°)(w^)  (Fig.  A. 3-3).  The  com¬ 
bined  data  transformation  in  Fq .  A. 3-15  can  be  implemented  in 
order  N  log2  N  operations  since  multiplications  by  H  and  W  are 
both  order  N. 
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R-4*24« 


Figure  A. 3-3  Product  Window  Definition 


Corresponding  to  the  frequency  domain  data  vector  z' 
(consisting  of  Fourier  coefficients)  is  the  auto-covariance 
matrix 

C'zz  =  E{z  '  (z '  )T}  =  (H  F  W)  Czz  (H  F  W)T  (A.  3-16) 

Similarly,  introducing  the  transformed  estimates 

x'  =  (H  F  W)  x  (A. 3-17) 


leads  to  the  cross-covariance  matrix 

C’  =  E{x' (z ' )T}  =  (H  F  W)  C  (H  F  W)T  (A. 3-18) 

J\Z 

A  straightforward  row  and  column  implementation  of  Eqs .  A. 3-16 

2 

and  A. 3-18  would  result  in  an  order  N  log2  N  algorithm.  This  is 
avoided  by  replacing  C^z  and  Czz  by  their  banded  approximations, 
and  calculating  only  those  elements  within  the  retained  bands. 

If  the  number  of  bands  in  the  one-dimensional  window  is  this 


A-47 


THE  ANALYTIC  SCIENCES  CORPORATION 


reduces  the  workload  to  the  order  of  mg  N  log2  N>  where  for 
typical  applications  m0  is  small  and  independent  of  N.  The 

D 

scheme  used  to  compute  C’  and  C'  bandwise  is  derived  in 

K  xz  zz 

Ref.  16  and  is  an  essential  part  of  the  GEOFAST  algorithm. 

The  preceding  transformations  constitute  the  first 

stage  of  the  calculation  (Fig.  A. 3-2).  Note  that  storage  re- 

2 

quirements  for  the  covariances  are  of  order  rag  N.  The  second 
stage  consists  of  the  solution  of  the  estimation  equations  in 
the  frequency  domain.  The  estimation  equations  in  the  trans¬ 
formed  variables  become  simply 


=  C’  (C’  ) 
xz  zz 


(A. 3-19) 


coupled  with  the  inverse  of  Eq.  A. 3-17 


x  =  (H  F  W)"1  x' 


(A. 3-20) 


It  is  easily  verified  that  the  solution  of  Eqs.  A. 3-19  and 
A. 3-20  is  identical  to  the  solution  of  Eq .  A. 3-1.  As  before, 
the  equations  are  solved  in  two  steps 


(c;_  r1  z 

zz  — 


(A. 3-21) 


x’  =  C’  u’ 
—  xz  - 


(A. 3-22) 


Since  C'  is  banded,  the  multiplication  in  Eq .  A. 3-22  requires 

y*z  2 

only  on  the  order  of  mg  N  operations. 


However,  analysis  of  the  solution  of  Eq .  A. 3-21  is 
more  complicated.  It  is  shown  in  Ref.  16  that  the  block 
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Toeclitz  structure  of  C zz  in  Eq.  A. 3-8  gives  rise  to  a  corres¬ 
ponding  block-banded  structure  for  C'  in  which  the  blocks  are 

zz 

also  banded  (Fig.  A. 3-4).  There  is  an  essential  difference 
between  this  structure  and  the  simple  band  structure  of  a 
diagonal  block  which  is  the  one-dimensional  analog  of  C^z .  In 
the  one-dimensional  problem  the  symmetric  band  matrix  C'  can 

9  Z  Z 

be  inverted  in  order  nig  n  operations  by  the  standard  Cholesky 
technique  (Refs.  29,  30),  where  n  is  the  matrix  dimension  and 
mg  the  number  of  superdiagonal  bands.  When  the  Cholesky  tech¬ 
nique  is  applied  to  the  block-banded  structure  of  the  two- 
dimensional  problem,  the  internal  banding  within  the  blocks  is 
lost  during  the  calculation.  This  difficulty  is  known  as  "fill- 
in."  The  result  is  an  effective  bandwidth  of  mBn2 ,  where  mB 
is  now  the  number  of  blockbands  of  size  n2  by  n2 ,  and  a  workload 
of  order  (mBn2 ) ^  (n-^n2 ) .  Assuming  n^  =  n2  gives  an  order  of 

effort  proportional  to  n^  =  Nz  rather  than  N. 
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Figure  A. 3-4 
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If  it  is  desired  to  reduce  the  computing  effort  further, 
an  order  N  algorithm  for  the  solution  of  Eq.  A. 3-21  can  be 
developed  by  using  an  iterative  method  based  on  a  special  class 
of  covariance  matrices.  This  special  type  of  covariance  struc¬ 
ture  is  termed  separable  and  is  characterized  by  a. covariance 
function  in  Eq.  A. 3-14  of  the  form 

f zz(® 2”®  1  *  $2~^1  ^  =  ^1^®2”^1^  ^2^2"^1  ^  (A. 3-23) 

This  class  of  functions  is  not  general  enough  to  include  the 
statistical  models  required  for  gravimetric  applications,  but 
it  can  serve  as  an  approximation  for  use  in  an  iterative  method. 
It  is  shown  in  Ref.  16  that  for  a  separable  covariance  matrix 
Eq.  A. 3-21  "factors"  into  two  one-dimensional  matrix  equations 
of  the  form 


V'  =  (C^)'1(Z,)T 


(A. 3-24) 


U’  =  (C|)'1(V’)T  (A. 3-25) 

where  correspond  to  f  ^ and  are  square  band  matrices 

of  dimension  n^ .  ^2  *  The  matrices  Z'  and  U'  are  both  of  dimen¬ 
sion  n-^  by  n2 ,  and  are  equivalent  to  the  vectors  z'  and  u'  when 
their  elements  are  listed  in  row-by-row  order.  If  the  bandwidth 
of  the  matrix  C2 '  is  mg,  then  Eq.  A. 3-24  can  be  solved  by  the 
Cholesky  method  in  order  mgn2(nig+n^)  operations.  Assuming 
nl  =  n2  ’  t^ie  toCa^  workload  for  both  equations  is  thus  of  order 

mgN  in  the  separable  case. 

The  iterative  solution  of  Eq.  A. 3-21  in  the  nonseparable 
case  proceeds  as  follows.  Choose  the  matrix  D'  as  a  separable 
approximation  to  C'  (which  may  be  done  in  a  natural  way)  and 
define  the  matrix  E'  by 


A-50 


THE  ANALYTIC  SCIENCES  CORPORATION 


C'  =  D*  +  E r  (A. 3-26) 

zz 

The  matrix  E'  will  be  block-banded  since  both  C'  and  D'  are. 

zz 

The  estimation  equation  can  now  be  written  implicitly 

u'  =  (D* )_1  [z'  -  E'u' ]  (A. 3-27) 


or  in  recursive  form 

u'(k+l)  =  (D')"1  [z'  -  E'u'(k)]  (A. 3-28) 

where  k  indicates  the  iteration  number  and  u'(0)  is  any  conven¬ 
ient  initial  solution. 

Since  both  the  inversion  and  multiplication  in  Eq.  A. 3-28 
require  on  the  order  of  N  operations,  so  does  one  step  of  the 
iteration.  If  the  number  of  steps,  mg ,  necessary  for  convergence 
of  the  iteration  is  independent  of  N  and  small  enough,  the  total 
workload  to  solve  Eq.  A.  3-21  is  of  order  N.  It  can  be  shown 
that  the  rate  of  convergence  depends  on  the  norm  of  the  matrix 
(D1)  ^E ' ,  which  will  be  small  if  the  separable  approximation  is 
a  good  one.  Furthermore,  a  simple  modification  to  the  iteration 
guarantees  convergence  for  any  positive  definite  matrix  D'.  The 
convergence  properties  of  the  iterative  technique,  Eq.  A. 3-28, 
and  the  appropriate  choice  of  the  approximating  matrix  D'  are 
discussed  in  Ref.  16. 

The  solution  of  the  estimation  equations  in  the  fre¬ 
quency  domain  completes  the  second  stage  of  the  GEOFAST  algorithm 
(Fig.  A. 3-2).  The  third  stage  consists  of  the  inverse  trans¬ 
formation  back  to  the  space  domain  (Eq.  A. 3-20).  Since  W  is  a 
diagonal  matrix,  its  inverse  is  an  order  N  operation.  The 
complex  DFT  matrix  F  can  bt  inverted  in  order  N  log2  N  operations 
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by  means  of  the  inverse  FFT,  while  the  inverse  of  the  matrix  H 
can  be  obtained  analytically  and  has  the  same  sparse  structure 
as  H  itself.  Thus,  the  final  stage  has  a  workload  dominated 
by  N  log2  N. 


The  computer  time  and  storage  requirements  for  the 

GEOFAST  algorithm  are  summarized  in  Table  A. 3-1.  The  numbers 

given  are  asymptotic  estimates,  valid  for  large  N,  and  do  not 

include  multiplicative  factors  independent  of  N  except  for  m^ 

and  mg .  For  repeated  application  of  the  algorithm  with  the 

same  covariance  model  the  matrices  C'  and  C'  can  be  stored, 

xz  zz 

reducing  the  time  requirement  to  the  order  of  N  log2  N  for 
each  application  to  a  new  data  set. 


TABLE  A. 3-1 


COMPUTER  TIME  AND  STORAGE  REQUIREMENTS 


ALGORITHM 

ASYMPTOTIC  ORDER 

OF  MAGNITUDE 

STAGE 

— 

TIME 

STORAGE 

1: 

— 

Input 

Transforms 

2 

mB  N  log2 

N 

mB2» 

2: 

Estimation 

Equations 

2xr 

msmB  N 

«,b2n 

3: 

Output 

Transform 

N  log2  N 

N 

Full  Solution 

2 

nig  N  log2 

N 

mB2N 

N  =  Total  number  of  data  points 
mg  =  Bandwidth  in  one  dimension 

mg  =  Number  of  steps  for  convergence 
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A. 4  OTHER  FOURIER  TRANSFORM  METHODS 

When  the  flat-earth  approximation  is  valid  (Ref.  31), 
methods  based  on  the  two-dimensional  Fourier  transform  lead  to 
efficient  algorithms  for  a  variety  of  geodetic  computations. 

The  flat-earth  versions  of  the  Stokes  integral  for  geoid  undula¬ 
tion,  the  Vening  Meinesz  integral  for  deflections  of  the  verti¬ 
cal,  and  the  Poisson  integral  for  upward  continuation  are 
naturally  implemented  by  Fourier  transforms,  because  each  has 
the  form  of  a  two-dimensional  linear  convolution.  In  addition, 
the  flat-earth  terrain  correction  can  be  approximated  as  a  sum 
of  linear  convolution  integrals. 

This  appendix  considers  the  basic  theory  of  Fourier 
transform  methods.  Section  A. 4.1  states  the  conventions  and 
notation  used  for  two-dimensional  Fourier  transforms  and  convo¬ 
lutions.  The  Fourier  transform  implementations  of  the  Stokes 
and  Vening  Meinesz  integrals  are  outlined  in  Section  A. 4. 2. 
(Upward  continuation,  long  known  to  be  very  straightforward  in 
the  transform  domain,  is  omitted  from  further  discussion.) 

All  Fourier  transform  methods  have  several  elements 
in  common.  In  each  case,  the  basic  strategy  is  as  follows: 

•  Apply  appropriate  weighting  and/or  en¬ 
closing  windows  to  the  data 

•  Take  the  numerical  Fourier  transform  of 
a  set  of  measurements 

•  Multiply  by  the  appropriate  transfer 
function  (almost  always  computed  directly 
from  an  analytic  expression) 

•  Finish  by  returning  to  the  space  domain 
with  a  numerical  inverse  transform. 
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The  methods  are  fast,  because  they  take  advantage  of  the  speed 
of  the  Fast  Fourier  Transform  (FFT)  algorithm,  and  simultane¬ 
ously  compute  all  output  values.  They  require  a  regular  grid 
of  input  data,  and  return  output  values  on  the  same  geographical 
grid  as  the  input  data.  (Values  at  intermediate  or  arbitrary 
points  must  be  interpolated.)  As  with  most  Fourier  transform 
applications,  care  is  required  to  apply  reasonable  windowing 
to  the  data  (to  avoid  undesirable  edge  effects)  and  to  avoid 
spatial  aliasing  in  the  input  data. 

A. 4.1  Fourier  Transforms  and  Convolution  Integrals 

In  the  following  two  sections  Bracewell's  (Ref.  32, 
pp.  241-243)  conventions  for  the  two-dimensional  Fourier  trans¬ 
form  and  two-dimensional  convolution  integral  are  used.  The 
Fourier  transform  of  a  two-dimensional  function  f(x,y)  is  de¬ 
fined  by 

F(u,v)  =  JJ  f(x,y)e  (ux+vy ^xdy  (A. 4-1) 

where  x  and  y  are  Cartesian  coordinates  in  the  plane,  u  and  v 
are  Cartesian  coordinates  in  the  transform  domain,  and  i  =  J-l. 
The  integration  extends  over  the  entire  plane;  the  limits  are 
omitted  for  convenience.  The  inverse  transform  is  given  by 

f(x,y)  =  JJ  F(u  ,v)e^n  ^UX+V^dudv  (A.  4- 2) 

Here  again,  and  throughout  the  remainder  of  the  section,  inte¬ 
gration  limits  of  -®  and  +®  are  omitted  for  convenience. 

The  convolution  of  a  pair  of  two-dimensional  func¬ 
tionals  f(x,y)  and  g(x,y)  is  expressed  by  the  integral 
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f*g 


x,y 


f(x'  ,y '  )g(x-x'  ,y-y  '  )dx'dy  ' 


(A. 4-3) 


where  the  primes  denote  auxiliary  variables  of  integration. 
According  to  the  Convolution  Theorem  the  Fourier  transform  of 
a  convolution  is  simply  the  product  F(u , v)G(u , v) ,  where  G  de¬ 
notes  the  Fourier  transform  of  g.  Applying  Eq .  A. 4-2  yields 


f*g 


x,y 


ff  F<u,v)G<u,v)ei2*<“+v5')dudv 


(A. 4-4) 


which  is  the  principal  basis  for  the  methods  described  in  the 
remainder  of  this  section. 

A. 4. 2  Stokes  and  Vening  Meinesz  Integrals 

Derivation  of  the  Fourier  transform  versions  of  the 
flat-earth  Stokes  and  Vening  Meinesz  integrals  is  straightfor¬ 
ward,  because  there  is  no  need  for  approximations  to  cast  them 
in  the  form  of  linear  convolution  integrals;  the  exact  formulas 
already  are  in  such  form.  Thus,  as  outlined  in  Ref.  31  and 
implemented  by  Sideris  (in  Ref.  21)  and  at  TASC  (this  report, 
p.  A-57),  it  is  natural  to  apply  Fourier  transform  (hence  FFT) 
methods  to  the  computation  of  geoid  undulation  and  deflections 
of  the  vertical  from  gravity  anomalies. 

The  flat-earth  Stokes  integral  for  undulation  of  the 

geoid  is 

N(xp,yp)  =  ff  Ag-^y)  dxdy  (A. 4-5) 
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where  g  is  normal  gravity,  Ag  is  the  gravity  anomaly  function, 

and  rQ  =  [  (xp-x)  +(yp-y)'i]'5.  The  north  (meridional)  and  east 
(prime  vertical)  components  of  the  deflection  of  the  vertical, 
i  and  r|  ,  respectively,  are  given  by  the  flat-earth  Vening 
Meinesz  integral 


4(xp,y p)l  1 

(  =  2Si; 

n (Xp,yp)  I 


I  (yp-y)/rQ  I 

Ag(x,y)  1  (  dxdy 


(xp-x)/r0‘ 


(A. 4-6) 


where  the  y-axis  is  directed  north  and  the  x-axis  is  directed 
east.  (Note  that  |  =  -3N/3y  and  n  =  -3N/3x. ) 


The  Fourier  transforms  of  the  kernel  functions  in  Eqs .  A.  4-7 
and  A. 4-8  are  straightforward  closed-form  expressions.  The  FFT 
may  be  applied  to  gridded  gravity  anomaly  data  to  estimate  the 
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Fourier  transform  AG(u,v).  Thus,  one  may  evaluate  Eqs .  A. 4-7 
and  A. 4-8  efficiently,  at  all  points  on  a  regular  grid,  by 
multiplication  in  the  transform  domain  followed  by  an  inverse 
FFT ,  as  in  Eq.  A. 4-4.  The  relevant  Fourier  transform  pairs 
are  listed  in  Table  A. 4-1. 

Some  of  the  practical  aspects  of  using  geodetic  algo¬ 
rithms  based  on  FFT  implementation  of  the  convolution  Eqs.  A. 4- 7 
and  A. 4-8  have  already  been  mentioned.  In  the  late  1970s, 

TASC  implemented  computer  programs  using  a  mixed  radix  FFT  for 
performing  the  flat-earth  Stokes  and  Vening  Meinesz  computations. 
Work  with  this  software,  using  marine  gravity  and  altimetry 
data,  as  well  as  terrestrial  observations,  clarified  some  of 
the  advantages  and  disadvantages  of  the  Fourier  transform  ap¬ 
proach.  Unlike  collocation,  this  approach  is  deterministic, 
not  statistical.  An  advantage  relative  to  collocation  is  that 


TABLE  A. 4-1 

FOURIER  TRANSFORM  PAIRS 


.  FUNCTION 

TRANSFORM 

Ag(x,y) 

AG(u,v) 

2  2  ^ 
l/(xz+y^ ) 

2  2 

l/(u^+v  ) 

2  2  3/2 

Y/U  ) 

2  2  ^ 

-2n iv/( u^+v  ) 

2  2 
x/ (x^+y^ ) 

o  2  ^ 

- 2n iu/Cu^+v  ) 

/N(x ,y)\ 

lg 

AG( u , v )/( u2+v2  ) 

\  £<x,y) > 

\n(x,y)/ 

{  2  2 

i  |  v/(u^+v^>  ) 

AG(u,v)  J  ( 

•  °  }  2  2  M 

(u/(uz+vz)  ) 
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there  is  no  need  to  assume  covariance  models,  yet  the  corres¬ 
ponding  disadvantage  is  that  the  variances  and  covariances  of 
the  estimation  errors  are  not  available.  A  characteristic  of 
the  FFT  approaches  is  that  data  and  estimates  must  be  on  regular 
grids.  As  discussed  above  in  Section  A. 3.1  of  this  appendix, 
the  grid  requirement  is  consistent  with  practical  applications; 
the  computational  effort  of  gridding  the  data  is  more  than 
offset  by  the  much  improved  efficiency  of  these  methods. 

Other  important  aspects  of  the  FFT  approaches  that 
must  be  considered  are  the  need  to  control  spectral  leakage 
and  aliasing  by  proper  windowing  of  the  data  region.  TASC 
found  that  leakage  problems  could  be  controlled  by  first  remov¬ 
ing  a  low-order  spherical  harmonic  reference  field  from  the 
data,  then  adding  back  the  appropriate  reference  values  after 
applying  the  flat-earth  Stokes  or  Vening  Meinesz  algorithms  to 
the  reduced  data.  (This  approach  is  also  described  in  Ref.  21.) 
Edge  effects  have  been  successfully  controlled  in  work  at  TASC 
by  adding  a  continuous  border  region  to  the  data,  formed  by 
two-dimensional  reflection  of  the  outer  edge  of  the  data.  The 
frequency  domain  algorithms  are  then  applied  to  the  FFT  of  the 
extended  region. 

When  implemented  with  the  practical  modifications 
described  above,  the  FFT-based  flat-earth  Stokes  and  Vening 
Meinesz  algorithms  have  proven  to  be  very  fast,  accurate,  and 
easy  to  implement.  They  deserve  more  widespread  use. 


A- 58 


THE  ANALYTIC  SCIENCES  CORPORATION 


APPENDIX  B 

FFT  TERRAIN  CORRECTION  ALGORITHM:  DETAILS  OF 
NEW  DERIVATION 


This  appendix  describes  in  detail  the  new  derivation 
of  the  FFT  terrain  correction  algorithm  introduced  in  Chap¬ 
ter  5  and  Ref.  35.  As  noted  in  Chapter  5,  the  flat-earth  ter¬ 
rain  correction,  a  nonlinear  function  of  terrain  elevation 

h(x,y),  may  be  expressed  approximately  in  terms  of  linear  con- 

2 

volution  integrals  involving  h  and  h  .  The  appendix  begins 
with  a  review  of  the  previous  derivations  of  such  approximate 
forms  (Refs.  33  and  34).  It  is  followed  by  the  new  deriva¬ 
tion,  showing  how  the  approximation  is  more  rigorous  than 
claimed  by  its  proponents.  Finally,  two  analytic  examples 
illustrate  the  effects  of  various  forms  of  the  approximation. 

The  flat-earth  terrain  correction,  evaluated  at  the 
point  (xp ,yp ) ,  is  given  by  (e.g.,  Eq.  5-1;  Ref.  33,  Eq.  11; 
Ref.  34,  Eq.  7.16) 


c(xp,yp)  =  Gp 


If  dxdy  / 


b(x,y) 


U-hp) 


dz 


[(xp-x)2+(yp-y)2+(hp-z)2]3/2 


(B-l) 

where  G  is  the  gravitational  constant,  p  is  the  density  assumed 
for  the  terrain,  and  hp  =  h(xp,yp). 
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Now,  expanding  the  denominator, 


[ (xp-x)2+(yp-y)2+(hp-z)2] 


3/2 


2^  ^  ^  "  hp ) 


r3  2  r5 
o  o 


(B-2) 


2  2  is 

where  rQ  =  [ (xp-x)  +  (yp-y)  l"5.  Thus,  it  is  reasonable  to 
make  the  approximation 


1  _  1 


(B-3) 


2  2 

assuming  (z-hp)  /rQ  <<  1,  and  to  integrate  over  z  to  obtain 


c(xp  yp)  = 


$?Gp  If 


( h(x ,y ) -hp ] ' 


[  (xp-x)  +  (  yp-y )' 


3/2 


dxdy  (B-4) 


In  both  previous  derivations  (Ref.  33,  Eq .  14;  Ref.  34,  Eq.  7.17), 

this  is  an  important  intermediate  result.  Note  that  the  quantity 
2  2 

(z-hp)  /rQ  is  generally  small  for  large  ro  while  for  very  small 
rQ  it  is  bounded  by  the  square  of  the  topographic  slope  S: 

(z-hp)2/r2  <  (h(x,y)-hp)2/r2  <  S2  (B-5) 

Equation  B-4  may  be  expressed  as  a  sum  of  two  convo¬ 
lutions  plus  a  constant  term  by  expanding  the  numerator: 


c(xp,yp)  =  ~  Gp 

(h2*f) 

1  -  2h  (h*f) 

1  P 

+  hp2  if  f (x,y)dxdy 

1  xp,yp 

xp,yp  ^ 

B-2 
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2  2  -3/2 

where  f(x,y)  =  (x  +y  )  '  Equation  B-6  suggests  computing 

the  terrain  correction  by  applying  Eq .  A. 4-4  to  the  Fourier 

2 

transforms  (Eq.  A. 4-1)  of  f  and  h,  and  f  and  h  .  There  is  a 

2  2-3 /2 

problem,  however:  the  function  (x  +y  )  '  is  not  integrable 
over  the  plane  and  its  Fourier  transform  does  not  exist.  Two 
remedies  have  been  suggested. 

Reference  33  follows  the  approach  of  evaluating  f 
numerically  on  a  regular  x-y  grid,  replacing  the  singularity 
at  x=y=0  by  zero,  then  computing  a  numerical  Fourier  transform 
(using  the  FFT).  In  addition,  the  nonconvergent  integral 
appearing  in  the  third  term  in  Eq .  B-6  was  replaced  by  the 
value  of  the  numerical  transform  at  u=v=0 .  Obviously,  this  is 
an  ad  hoc  remedy,  and  it  is  difficult  to  evaluate  its  effect, 
since  a  sum  of  infinite  quantities  is  replaced  by  a  sum  of 
finite  quantities. 

The  second  remedy  was  suggested  in  Ref.  34  and  also 
mentioned  in  Ref.  33.  It  consists  of  an  analytic  modification 
of  f  such  that 

f(x,y)  =  (x2+y2+a2)'3/2  (B-7) 

I 

where  a  is  a  constant  that  is  chosen  empirically.  The  modified 
f  has  the  Fourier  transform 

F(u,v)  =  2na  ^e  27ta^U  +v  (B-8) 


and 


f(x,y) 


dxdy  = 


F  ( 0 , 0  ) 


2na 


-1 


( B-9 ) 
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Although  it  is  more  straightforward  to  implement,  avoiding  a 
numerical  transform  of  the  modified  f,  this  remedy  appears  as 
ad  hoc  as  the  first;  also  there  is  an  arbitrary  constant  to  be 
chosen . 


However,  Eq.  B-7  combined  with  Eq.  B-6  is  a  better 
approximation  to  Eq .  B-l  than  is  Eq.  B-4,  given  a  suitable 
choice  of  the  constant  a.  The  FFT  algorithm  for  terrain  ef¬ 
fect  computation  is  actually  more  rigorous  than  its  proponents 
have  claimed. 

To  understand  why  this  is  so,  consider  the  following 
expressions : 


h(x,y) 


cUp,yp)  ^  Gp 


ff  dxdy  / 


(z-hp)dz 


[(xp-x)2-Kyp-y)2+a2] 


3/2 


(B-lOa) 


=  >// 


[h(x,y)-hp]‘ 


[(Xp-x)  +(yp-y)  +a" 


3/2 


dxdy 


(B-lOb) 


Equation  B-lOb  is  equivalent  to  the  combination  of  Eqs.  B-6 
and  B-7;  Eq.  B-lOa  is  identical  to  Eq .  B-l  except  that  in  the 
denominator  (hp-z)  is  replaced  by  a  .  Thus,  the  combination 
of  Eqs.  B-6  and  B-7  follows  directly  from  the  exact  (flat-earth) 

terrain  correction  (Eq.  B-l),  given  the  single  approximation 

2  2  2  -  - 

(hp-z)  =  a  .  The  meaning  of  a  is  now  clear;  it  represents 

v  2 

an  average  value  of  (hp-z)  and  its  optimal  value  always  ex¬ 
ceeds  zero  because  (hp-z)^  is  always  positive. 
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Thus,  the  two  stages  of  approximation  applied  in 

Refs.  33  and  34  are  unnecessary  and  partially  cancel  one 

3 

another.  Truncation  of  the  series  expansion  of  1/r  leads  to 
an  overestimate  of  the  terrain  correction  which  is  compensated 
by  the  fact  that  (x^+y^+a^)  underestimates  (x^+y^)  . 

Examples  with  exact  analytic  solutions  help  to  illus¬ 
trate  and  evaluate  the  various  approximations  to  the  flat-earth 
terrain  correction.  Two  such  examples  are  given  here.  The 
first  involves  the  conical  mountain  used  in  Ref.  34  (pp.  46-47). 
The  second  is  based  on  an  "inner  zone"  of  constant  slope. 
(Note  that  the  approximations  in  the  integrands  of  Eqs .  B-4 
and  B-lOb  are  both  most  severe  when  rQ  is  small.) 

The  conical  mountain,  with  height  H  and  slope  angle 
6,  is  illustrated  schematically  in  Fig.  B-l.  The  terrain  cor¬ 
rection  at  the  summit  point  P  is  straightforward  to  evaluate. 
Placing  the  origin  such  that  xp=yp=0 ,  the  height  function  is 

(  H  -  r  tan  0  r  <  H/tan  0 

h(x,y)  =  h(r  )  =  °  (B-ll) 

(0  ro  -  H/tan  0 

f~~2  2 

where  rQ  -  Jx  +y  .  Thus,  hp  =  H.  Substitution  into  Eq .  B-l 
gives  the  exact  result 


;t  =  2nGp 


/  rodro  / 


h(ro> 


(z-H)dz 


l rQ  + (H-z ) ^ ] 


=  2nGpH  sin  0 


(B-12) 
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P 


Figure  B-l  Conical  Mountain 


(as  shown  in  Ref.  34,  Eq.  7.24).  On  the  other  hand,  substitu¬ 
tion  into  Eq.  B-lOb  yields  the  approximate  terrain  correction 


capprox<a)  =  *<30 


/ 


“  [h(ro)-H]2 


[ro2+a2] 


3/2 


r  dr 
o  o 


2  2  2 

=  2nGp  tan  0  [ (H  +a  tan  0)  -a  tan  0] 


(B-13 ) 

which  is  a  function  of  the  parameter  a.  Substitution  of  h(x,y) 
into  Eq.  B-4  is  equivalent  to  setting  a=0 , 

c  (0)  =  2nGpH  tan  0  (B-14) 

approx 

in  agreement  with  Eq.  7.25,  Ref.  34. 

2 

As  pointed  out  earlier,  the  optimal  value  of  a  ex¬ 
ceeds  zero.  In  fact,  for  the  single  point  P,  a  may  be  chosen 
such  that  capprox(a)  =  cexact.  By  solving  the  equation 
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¥ 


tan  0  1 (H2+a2tan20 )  -a  tan  0]  =  H  sin  0  (B-15) 

for  a,  one  finds  a  =  H  (sin  0)/2.  Of  course,  setting  a  = 

H  (sin  0  )/2  in  Eq .  B-lOb  probably  would  not  lead  to  exact  re¬ 
sults  for  points  on  the  flank  of  the  cone,  or  on  the  surround¬ 
ing  plain,  but  the  point  is  that  a  range  of  nonzero  values  of 
a  does  better  than  a=0.  Furthermore,  for  the  summit  point  P, 
one  can  derive  an  upper  bound  of  that  range  by  setting  the 
underestimate  of  the  terrain  correction  incurred  by  using 

c=.nr.T-r.v(a)  equal  to  the  overestimate  incurred  using  c  (0): 

approx  &  approx 

cexact  capprox^a^  “  capprox^^  "  cexact 

H  sin  0  -  tan  0  [ (H2+a2tan20 )  -a  tan  0 ]  =  H  tan  0  -  H  sin  0 


2H  cos20  (1  -  cos  0) 
sin  0  (2  cos  0-1) 


(B-16) 


Table  B-l  presents  the  computed  values  of  these  results  for 
several  values  of  0  . 


I 

I 

I 


TABLE  B-l 


CONICAL  MOUNTAIN  EXAMPLE 


0 

BEST  a 

UPPERBOUND  ON  a 

Small  0  (in  radians) 

*50  H 

0H 

15  deg 

0.13H 

0.26H 

30  deg 

0.25H 

0.55H 

45  deg 

0.35H 

H 
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The  table  indicates,  for  instance,  that  for  0-30  deg 
c  =  c  _  if  a  =  0.25H  and  any  value  of  a  between  0 

and  0.55H  yields  a  better  approximation  than  does  a  =  0. 

Another  straightforward  example  can  be  developed  by 

considering  the  inner  zone  in  which  rQ  is  small  enough  to  allow 

the  approximation  of  the  terrain  correction  by  its  first-order 

Taylor  expansion.  This  is  the  region  where  the  approximations 

to  1/r  have  the  greatest  effect.  Consider  the  case  where  0  5 

r  <  R  and 
o  — 


xP  =  yp  =  0 

h(x,y)  =  h(rQ,<j>)  =  hp  +  S  rQcos  $ 


(B-17 ) 


where  x  =  r  cos  0  and  y  =  r  sin  0  .  This  defines  a  circular 
o  o 

zone  of  constant  slope  S  with  its  gradient  directed  (arbitra¬ 
rily)  in  the  x-direction.  From  Eq.  B-l,  the  exact  terrain 
correction  due  to  this  zone  is 


'exact 


=  GpR 


2n 

/ 


[1-(1+S2cos20  )  2)  d0 


=  Jgps2r  <1  -  xe  s2  +  • • • ) 


(B-18) 
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Equation  B-10  yields 


C  (a) 

approx 


-h 


-2a] 


(B-19) 


and  setting  a  =  0  or  applying  Eq.  B-4  yields 


C  (0) 

approx 


( B-20 ) 


Clearly,  the  approximate  solution  can  be  closely  matched  to  the 

9  2 

exact  by  choosing  a  =  R  while  for  a  slope  corresponding  to 

30  degrees,  c  (0)  includes  an  error  of  almost  20%. 

app rox 

The  analytic  examples  presented  here  complement  the 
numerical  example  given  in  Chapter  5,  and  along  with  details 
of  the  new  derivation  reaffirm  the  conclusions  that  the  FFT 
terrain  correction  algorithm  is  more  rigorous  than  claimed  by 
its  proponents  and  that  the  optimal  value  of  the  parameter  a 
always  exceeds  zero. 


THE  ANALYTIC  SCIENCES  CORPORATION 


APPENDIX  C 

ANALYSIS  OF  TERRAIN  EFFECTS  USING 
WHITE  SANDS  DATA 


C . 1  INTRODUCTION 

In  a  region  of  high  relief  such  as  the  White  Sands 
study  area,  it  is  important  to  include  the  effects  of  terrain 
in  the  prediction  of  gravity  anomaly  or  deflections  of  the 
vertical  from  gravity  anomaly  observations.  The  research  de¬ 
scribed  in  this  appendix  addresses  the  quantitative  determina¬ 
tion  of  the  importance  of  terrain  effects,  including  dependence 
on  wavelength,  variation  with  terrain  type,  and  the  orientation 
of  prominent  terrain  and  gravity  features.  This  complements 
the  work  of  the  IAG  Study  group  (Refs.  20,  21)  on  algorithms 
for  prediction  of  gravity  disturbance  and  deflection  of  the 
vertical . 


The  investigations  documented  in  this  appendix  are 
summarized  for  quick  overview  in  Fig.  C.l-1.  They  include  two 
distinct  studies.  One  of  these,  as  outlined  in  the  upper  half 
of  Fig.  C.l-1,  makes  use  of  spectral  techniques  to  examine  and 
compare  observed  free-air  gravity  anomalies  and  computed  topo¬ 
graphic  effects  for  two  types  of  terrain  --  mountain  and  pla¬ 
teau  --  and  for  two  profile  orientations  ---  east-west  and  north- 
south.  The  methodology  and  conclusions  of  this  study  are  the 
subject  of  Sections  C.2,  C.3,  and  C.4  of  this  appendix. 

The  second  study,  outlined  in  the  lower  half  of 
Fig.  C.l-1,  examines  how  well  deflections  of  the  vertical  can 
be  predicted  from  local  topographic  data,  supplemented  by  a 
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Figure  C.l-1  Overview  of  Terrain  Effects  Analysis 


long-wavelength  component  derived  from  a  global  geopotential 
model.  This  study  is  described  in  Section  C.5  of  the  appendix. 


The  first  step  in  the  spectral  studies,  described  in 
Section  C.2,  is  to  grid  and  display  the  free-air  gravity  anom¬ 
aly  data  over  a  3  deg  by  3  deg  portion  of  the  White  Sands  area, 
and  to  display  gridded  topographic  elevations  over  the  same 
3  deg  by  3  deg  area.  Six  profiles  of  gravity  and  topography 
are  selected  from  these  data  sets  and  gravity  effects  due  to 
topography  are  computed  for  each  profile  using  a  particular 
point-mass  algorithm  which  is  described  and  illustrated  in 
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Section  C.3.  Section  C.4  describes  the  estimation  of  power 
spectral  densities  and  the  coherence  between  free-air  gravity 
and  the  topographic  effect  from  selected  subsets  of  the  pro¬ 
files  (sorted  by  terrain  type  and  profile  orientation)  by  means 
of  state-space  modeling,  and  includes  illustration  and  analy¬ 
sis  of  the  results.  (This  satisfies  the  need  for  spectral 
analysis  indicated  by  K.P.  Schwarz  in  the  1983  Study  Group 
Report  (Ref.  20).) 

In  the  second  study,  the  point  mass  algorithm  of  Sec¬ 
tion  C.3  is  applied  to  predict  the  topographic  contributions 
to  the  deflections  of  the  vertical  at  each  of  the  astrogeodetic 
truth  data  locations.  Using  these  predictions,  bounds  are 
placed  on  the  significance  of  the  topographic  effect.  The 
astrogeodetic  data  are  too  sparse  to  support  the  use  of  grid- 
ding  followed  by  state-space  modeling  like  that  applied  to  the 
gravity  data. 


C.2  GRIDDED  GRAVITY  AND  TOPOGRAPHIC  DATA 

Gravity  data  were  supplied  in  the  form  of  free-air 
point  gravity  anomaly  measurements  unevenly  distributed  over  a 
5.5  by  5.5  deg  area  (Ref.  20).  A  3  by  3  deg  subset  of  this 
area  extending  from  32  to  35  deg  North  and  105  to  108  deg  West 
was  selected  for  this  study.  This  subset  was  selected  because 
it  contained  a  relatively  dense  concentration  (3855  points)  of 
gravity  observations  and  was  also  covered  by  digital  gridded 
topographic  data.  Figure  C.2-1  shows  a  map  of  the  locations 
of  the  3855  gravity  measurements  used  in  this  study. 

In  order  to  perform  spectral  analysis  of  free-air  point 
gravity  measurements,  it  is  necessary  to  grid  the  data.  A  121 
by  121  grid  (40  grid  intervals  per  degree)  is  computed  over 
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Figure  C.2-1  Locations  of  Free-Air  Gravity  Anomaly 
Points  in  the  Study  Area 


the  study  area  using  the  Barnes  analysis  scheme  (Ref.  42). 
The  gridding  algorithm  consists  of  a  computationally  simple, 
Gaussian  weighted-averaging  technique  which  assigns  a  weight 
to  a  datum  solely  as  a  function  of  distance  between  datum  and 
grid  point. 


Figure  C.2-2  is  a  color  shaded-relief  image  of  the  3 
by  3  deg  free-air  gravity  anomaly  grid.  The  magnitudes  of  the 
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gravity  values  range  from  approximately  -60  to  165  mgals,  and 
are  indicated  by  hue,  as  shown  on  the  scale  bar.  Shading  is 
included  to  enhance  short-wavelength  features  in  the  data. 
The  shading  is  computed  using  a  southeast-to-northwest  direc¬ 
tional  derivative. 

The  gridded  topographic  data  covering  the  3  by  3  deg 
study  area  are  shown  as  a  color,  shaded-relief  image  in 
Fig.  C.2-3.  The  sampling  interval  of  the  topographic  data  is 
120  samples  per  degree  (a  360  by  360  grid).  This  grid  sam¬ 
pling  is  three  times  as  dense  as  the  sampling  interval  of  the 
gridded  gravity  and  accounts  for  the  greater  sharpness  of  the 
topography  image.  The  scale  bar  provides  a  conversion  between 
hue  and  elevation  as  in  the  gravity  image.  The  shading  is 
based  on  a  directional  derivative  simulating  illumination  from 
the  southeast. 


C . 3  PROFILES  OF  ANOMALIES  AND  TOPOGRAPHIC  CORRECTIONS 

Three  types  of  profiles  are  discussed  in  this  section: 
free-air  gravity  anomaly  profiles,  topographic  effect  profiles, 
and  residual  anomaly  profiles.  The  free-air  gravity  anomaly 
profiles  are  taken  along  rows  or  columns  of  the  free-air  grav¬ 
ity  anomaly  grid.  Topographic  effect  profiles  consist  of  grav¬ 
ity  disturbance  values  computed  from  an  algorithm  which  is 
discussed  below.  The  residual  anomaly  profiles  result  from 
subtracting  topographic  effect  profiles  from  the  corresponding 
free-air  gravity  anomaly  profiles. 

An  algorithm  which  computes  the  vertical  component  of 
the  gravity  disturbance  from  topographic  data  was  used  to  com¬ 
pute  the  topographic  effect  profiles.  Figure  C.3-1  illustrates 
the  essential  elements  of  this  algorithm.  The  topographic 


C-6 


THE  ANALYTIC  SCIENCES  CORPORATION 


LOCATION  Of  ELEVATION  *  son 


Figure  C.3-1  Pictorial  Representation  of  Gravity 
Disturbance  Computation 


data  are  treated  as  a  grid  of  vertical  rectangular  prisms  ex¬ 
tending  upward  from  sea  level.  The  mass  of  each  prism  is  com¬ 
puted  by  multiplying  the  volume  times  the  density,  and  is  con¬ 
sidered  to  be  a  point  mass  located  at  the  center  of  each  prism. 
The  vertical  component  of  the  gravity  disturbance  at  P  (see 
Fig.  C.3-1)  is  then  computed  by  summing  the  point  masses  over 
a  square  window  as  follows: 


6gp  = 


w 

E 


(C.3-1) 


where 


<5gp  =  the  vertical  component  of  the  gravity 
disturbance  at  point  P 

w  =  window  width 


C-8 


THE  ANALYTIC  SCIENCES  CORPORATION 


mass  of  each  prism 

gravi rational  constant 

2  .  2  ,  2 

x  +  y  +  z 

The  disturbance  values  closely  approximate  the  complete  Bouguer 
correction  because  both  slab  and  terrain  effects  are  taken 
into  account  by  setting  sea  level  as  the  lower  boundary  of  the 
prisms . 


The  profile  sampling  density  for  the  topographic  effect 
equals  that  of  the  gridded  topography,  120  points  per  degree. 

In  order  to  match  the  point  locations  and  sampling  density  of 
the  free-air  gravity  anomaly  profiles  (40  points  per  degree), 
the  profile  data  series  are  averaged  and  subsampled  using  a 
three-point  moving  window. 

Once  the  free-air  gravity  anomaly  profiles  and  topo¬ 
graphic  effect  profiles  consist  of  identical  point  locations 
and  spacing,  it  becomes  possible  to  take  the  difference  between 
the  two,  forming  the  residual  anomaly  profiles. 

A  total  of  six  profile  locations  are  selected  to  cover 
a  variety  of  terrain  types.  Figure  C.3-2  shows  the  profile 
locations  superposed  on  a  contour  map  of  the  White  Sands  top¬ 
ography.  Three  profiles  run  east-west,  and  are  labeled  EW-1, 

-2  and  -3,  and  three  run  north-south,  and  are  labeled  NS-1, 

-2 ,  and  -3 . 

At  each  of  the  profile  locations  shown  in  Fig.  C.3-2, 
the  free-air  gravity  anomaly,  topographic  effect  and  residual 
anomaly  profiles  are  computed.  These  three  profile  types, 
along  with  terrain  elevation  profiles,  are  plotted  for  each 
location  and  are  shown  in  Figs.  C.3-3  through  C.3-8.  Subsets 
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Figure  C.3-3  Profile  EW-1 
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Figure  C.3-4  Profile  EW-2 
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Figure  C.3-7  Profile  NS-2 
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MOUTH  TOPOGRAPHIC  EFFECT  PROFILES  SOUTH 


TERRAIN  ELEVATION 


Figure  C.3-8  Profile  NS-3 
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of  mountainous-  and  plateau-type  terrains  are  also  indicated 
on  these  plots.  The  predominantly  negative  values  of  the  re¬ 
sidual  anomalies  are  consistent  with  regional  Bouguer  gravity 
for  the  southwestern  U.S.  (Ref.  43). 

The  next  section  of  this  appendix  presents  further 
discussion  and  interpretation  of  these  profiles  in  terms  of 
power  spectral  densities  and  spectral  coherences,  derived  using 
modern  state-space  estimation  techniques. 


C . 4  POWER  SPECTRAL  DENSITIES  OF  PROFILES 

Segments  of  data  representing  mountainous-  and  plateau- 
type  terrains  are  selected  from  the  six  profiles  shown  in 
Fig.  C.3-2.  The  actual  segments  used  are  indicated  on  the 
profile  plots  in  Figs.  C.3-3  through  C.3-8.  The  criterion 
used  to  classify  the  terrain  into  these  two  types  defines  ele¬ 
vations  above  1750  m  being  defined  as  mountainous,  and  eleva¬ 
tions  below  1500  m  as  plateau. 

The  purpose  of  the  spectral  analysis  of  these  profile 
segments  is  to  determine  spectral  differences  between  terrain 
types  and  to  determine  the  extent  to  which  directional  (east- 
west  versus  north-south)  anisotropy  exists.  In  order  to  make 
these  spectral  comparisons,  the  data  segments  are  pooled  into 
four  groups:  1)  east-west  plateau,  2)  north-south  plateau, 

3)  east-west  mountain,  and  4)  north-south  mountain.  Power 
spectra  for  each  of  the  four  groups  are  computed  by  state- 
space  modeling,  as  described  below. 

A  method  of  state-space  modeling  (Ref.  44)  is  used  to 
produce  the  power  spectral  density  (PSD)  estimates  discussed 
in  this  section.  This  method  consists  of  two  steps.  First,  a 
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model  is  developed  for  the  data  series  using  a  canonical-variates 
(CV)  technique.  The  CV  technique  produces  a  family  of  optimal 
state-space  models  by  solving  a  finite  number  of  linear  equations 
for  the  state-space  model  parameters.  In  the  second  step,  the 
power  spectrum  of  the  model  is  computed.  Once  the  model  param¬ 
eters  are  computed,  many  statistical  analyses  may  be  performed, 
including  an  analysis  of  coherence  between  the  data  series. 

This  state-space  modeling  method  is  advantageous  for  our  pur¬ 
poses  because  it  allows  multiple  tracks  of  data  to  be  pooled 
for  spectrum  estimates  that  are  limited  by  the  total  extent  of 
the  data,  rather  than  the  lengths  of  the  individual  track 
segments . 


Figures  C.4-1  through  C.4-4  show  power  spectral  den¬ 
sity  (PSD)  plots  of  free-air  gravity  anomaly,  topographic  ef¬ 
fect,  and  residual  anomaly  for  each  of  the  four  groups  of  data: 
east-west  plateau,  north-south  plateau,  east-west  mountain, 
and  north-south  mountain.  Coherence  between  the  topographic 
effect  and  free-air  gravity  anomaly  data  series  is  also  plotted 
for  each  of  the  four  data  groups. 

C.4.1  Discussion  of  Results 

As  noted  above,  the  selected  data  segments  representing 
mountainous-  and  plateau-type  terrains  in  east-west  and  north- 
south  profile  orientations  are  pooled  into  four  distinct  groups. 
The  various  PSDs  estimated  from  these  four  groups  are  discussed 
below. 


East-west  Plateau  -  The  data  segments  used  to  form 
this  pool  consist  of  two  segments  from  profile  EW-2  (Fig.  C.3-4) 
and  one  segment  from  profile  EW-3  (Fig.  C.3-5).  Figure  C.4-1 
shows  the  PSDs  for  the  free-air  gravity  anomaly,  topographic 
effect,  and  residual  anomaly.  The  PSD  plot  shows  that  at  all 
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Figure  C.4-3  East-West  Mountainous  Terrain 
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Figure  C.4-4  North-South  Mountainous  Terrain 
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frequencies,  the  topographic  effect  has  less  power  than  the 
free-air  gravity  anomaly  or  residual  anomaly.  A  local  peak  in 
power  at  40  km  wavelength  is  seen  in  all  three  PSDs.  This 
peak  may  be  due  to  the  large  north-south-trending  geologic 
features  crossing  the  White  Sands  study  area  (see  Figs.  C.2-2 
and  C.2-3)  . 

The  plot  of  coherence  between  free-air  gravity  anomaly 
and  topographic  effect  data  series  shows  that  there  is  low 
coherence  at  longer  wavelengths  (>60  km)  and  higher  coherence 
at  shorter  wavelengths,  particularly  in  the  region  of  the  local 
peak  of  power.  The  sudden  drop  in  coherence  near  the  folding 
frequency  is  due  to  the  finite  sample  spacing  and  is  not  of 
geophysical  significance. 

North-south  Plateau  -  The  data  segment  used  for  the 
north-south  plateau  group  is  from  profile  NS-2  (Fig.  C.3-7). 

The  plot  of  the  three  PSDs  (Fig.  C.4-2)  shows  that  the  topo¬ 
graphic  effect  has  very  low  power.  Subtraction  of  the  topo¬ 
graphic  effect  from  the  free-air  gravity  anomaly  does  not  sig¬ 
nificantly  change  the  amount  of  power  in  the  gravity  anomaly 
PSD;  terrain  plays  a  very  minor  role  in  this  case.  The  coher¬ 
ence  between  the  free-air  gravity  anomaly  and  topographic  effect 
PSD  is  extremely  low  at  all  frequencies. 

East-west  Mountain  -  The  data  segments  used  to  form 
this  group  are  from  two  segments  of  profile  EW-1  (Fig.  C.2-3) 
and  one  segment  of  profile  EW-3  (Fig.  C.3-5).  Examination  of 
the  PSD  plot  (Fig.  C.4-3)  shows  a  local  peak  in  power  similar 
to  that  observed  in  the  E-W  plateau  group.  Unlike  the  E-W 
plateau  group,  though,  the  peak  appears  only  in  the  topographic 
effect  and  free-air  gravity  anomaly  PSDs,  and  is  absent  in  the 
residual  anomaly  PSD.  Since  the  free-air  and  topographic  ef¬ 
fect  PSDs  track  each  other  very  well,  much  of  the  power  in  the 
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free-air  gravity  anomaly  data  is  due  to  terrain.  A  large  de¬ 
crease  in  the  amount  of  power  of  all  three  PSDs  at  short- 
wavelengths  is  also  observed. 

The  peak  in  coherence  corresponding  to  the  PSD  power 
peak  is  located  at  a  wavelength  of  about  40  km,  similar  to  the 
coherence  plot  in  Fig.  C.4-1.  The  coherence  between  free-air 
gravity  anomaly  and  topographic  effect  is  fairly  high  at  most 
frequencies;  therefore,  along  these  profiles  most  of  the  grav¬ 
ity  signal  reflects  topography  directly. 

North-south  Mountain  -  The  data  segments  used  to  form 
this  group  consist  of  one  segment  from  profile  NS-1  (Fig.  C.3-6), 
one  segment  from  NS-2  (Fig.  C.3-7),  and  one  segment  from  NS-3 
(Fig.  C.3-8).  Figure  C.4-4  shows  the  three  PSDs  computed  for 
this  group.  The  power  in  ail  three  PSDs  rolls  off  (decreases) 
with  a  decrease  in  wavelength.  The  coherence  begins  to  decrease 
sharply  in  the  vicinity  of  this  roll-off  of  power.  At  long 
wavelengths  (>100  km)  there  is  high  coherence;  at  shorter  wave¬ 
lengths,  the  coherence  drops  off  dramatically. 

C.4.2  Conclusions 

The  PSD  and  coherence  results  of  the  preceding  section 
are  presented  in  Table  C.4-1.  These  results  strongly  reflect 
the  anisotropy  in  the  terrain:  east-west  profiles  show  higher 
correlation  between  terrain  and  free-air  gravity  anomaly  than 
north-south  profiles.  This  is  because  the  predominant  struc¬ 
tural  features  of  the  terrain  (mountain  ridges  and  valleys) 
run  north-south.  The  correlation  is  strongest  at  wavelengths 
between  25  and  50  km. 

This  conclusion  reinforces  the  importance  of  removing 
terrain  effects  before  applying  collocation  techniques.  The 
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TABLE  C. 4-1 

COMPARISON  OF  MOUNTAIN  AND  PLATEAU  REGIONS  IN  NORTH- SOUTH 
AND  EAST-WEST  DIRECTIONS 


FREE-AIR 

GRAVITY  ANOMALY 

TOPOGRAPHIC  EFFECT 

_ 

RESIDUAL  ANOMALY 

DATA  TYPE 

RMS 

( mga 1 ) 

- 3dB  WAVELENGTH 
(km) 

— 

RMS 
(mgal ) 

-3dB  WAVELENGTH 
(  km ) 

RMS 
(ragal ) 

-3dB  WAVELENGTH 
(km) 

- —  ■ 

East -West  PLateau 

16.3 

23 

6.0 

17 

12.0 

19 

North-South  Plateau 

5.8 

65 

2.0 

100 

5.6 

50 

East -West 

Mount  a i n 

13.2 

25 

13.6 

25 

7.2 

20 

No  r  t  h  -  .south 

Mount  i  i n 

_ 

21.5 

200 

16.  A 

200 

12.9 

36 

SQUARED 

COEFFICIENT 

CORRECTION 

0.33 

0.05 

0 . 67 

0.62 


need  to  do  this  is  obvious  for  mountainous  areas;  it  was  pointed 
out  by  the  IAG  study  group.  The  present  results,  however, 
indicate  that  even  in  plateau  areas,  where  the  topographic 
contribution  to  free-air  gravity  anomaly  is  small,  there  is 
substantial  correlation  between  topography  and  gravity  anomaly 
for  the  east-west  tracks.  Hence,  for  such  areas,  it  will  be 
worthwhile  to  compensate  for  possible  anisotropy  by  applying 
terrain-remove-and-restore  techniques  before  carrying  out  col¬ 
location  or  other  estimation  procedures. 


C . 5  TERRAIN  EFFECTS  ON  DEFLECTIONS  OF  THE  VERTICAL 

A  control  data  set  containing  441  astrogeodetic  obser¬ 
vations  of  deflections  of  the  vertical  was  provided  in  the  IAG 
SSG  C.70  test  data  set,  version  la  (Ref.  20).  These  data  are 
located  primarily  in  the  central  region  of  the  3  deg  by  3  deg 
study  area  (32  deg  to  35  deg  North  and  105  deg  to  108  deg  West). 
Because  the  data  are  not  dense  enough  to  support  the  computation 
of  a  regular  grid,  the  spectral  analysis  method  used  in  the 
previous  section  is  not  feasible.  The  present  analysis  examines 
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the  prediction  of  the  irregularly  spaced  astrogeodetic  deflec¬ 
tions  as  the  sum  of  a  short-wavelength  component  (computed 
directly  from  the  terrain  height  data  using  a  point-mass  algo¬ 
rithm)  and  a  long-wavelength  component  derived  from  a  global 
geopotential  model.  An  improved  prediction  method  can  include 
contributions  from  measured  gravity  anomaly  data,  using  algo¬ 
rithms  such  as  Vening  Meinesz  integration  or  least-squares 
collocation  (Appendix  A;  Refs.  20  and  21),  although  here,  as 
in  Section  C.l,  the  terrain  contribution  is  emphasized. 

C.5.1  Discussion  of  Results 

Figure  C.5-1  shows  the  441  observed  deflections  of 
the  vertical  used  in  this  analysis.  In  this  figure,  the  arrows 
indicate  the  direction  and  magnitude  of  the  deflection,  with 
the  location  of  the  deflection  station  shown  at  the  tail  of 
the  plotted  vector.  It  is  important  to  note  that  the  east 
components  of  the  deflections  of  the  vertical  are  opposite  in 
sign  to  those  plotted  in  the  most  recent  IAG  study  group  report 
(Ref.  21,  Fig.  1.6,  p.  11).  Earlier  work  (Ref.  18,  Fig.  10, 
p.  7852)  showed  a  subset  of  the  observed  data  with  the  east 
components  of  the  deflections  of  the  vertical  in  agreement 
with  those  shown  in  Fig.  C.5-1.  Investigation  of  this  dis¬ 
crepancy  has  revealed  that  the  data  set  distributed  on  tape  is 
erroneous.  The  east  component  of  the  deflection  of  the  verti¬ 
cal  (n)  is  generally  much  larger  than  the  north  component  (£), 
largely  as  a  result  of  the  north-south  trending  topographic 
features  present  in  the  study  area. 

The  point-mass  algorithm  used  to  predict  gravity  dis¬ 
turbance  (see  Eq .  C.3-1)  is  modified  to  permit  computation  of 
the  north  and  east  components  of  the  deflections  of  the  ver¬ 
tical.  Using  the  gridded  topographic  data  described  in  Sec¬ 
tion  C.2,  the  point-mass  algorithm  is  applied  to  compute 
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LONGITUDE  (W) 

Figure  C.5-1  Observed  Deflections  of  the  Vertical 


deflection  values  at  each  of  the  441  observed  data  locations. 
The  resulting  terrain-predicted  deflections  of  the  vertical 
are  shown  in  Fig.  C.5-2.  These  values  represent  a  significant 
portion  of  the  short-wavelength  component  of  the  observed  de¬ 
flections  of  the  vertical. 

Figure  C.5-3  displays  deflections  of  the  vertical 
computed  from  Rapp's  1978  (Ref.  37)  set  of  geopotential  coeffi¬ 
cients  to  degree  and  order  180,  at  the  observed  data  locations. 


Figure  C.5-2  Terrain-Predicted  Deflections  of  the  Vertical 


expansion  of  the  coefficients  and  evaluate  the  deflection  com¬ 
ponents  on  a  regular  grid  in  colatitude  and  longitude  (Ref.  45) 
The  deflections  predicted  from  the  geopotential  model  are  fair¬ 
ly  constant  in  both  magnitude  and  direction;  they  approximate 
the  mean  of  the  observed  deflections. 


The  long-wavelength  contribution  predicted  by  the 
Rapp  1978  model  was  compared  with  two  other  prediction  models: 
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Figure  C.5-3  Deflections  of  the  Vertical  from 
Geopotential  Model 


•  The  long-wavelength  deflections  of  the 
vertical  computed  from  Rapp  1981  (Ref.  38) 
geopotential  model 

•  The  long-wavelength  components  obtained 
from  a  low  degree  polynomial  fit  to  the 
observed  data. 


The  Rapp  1981  model  produces  deflections  that  are  nearly  the 
same  as  those  computed  from  the  Rapp  1978  (Ref.  37)  coeffi¬ 
cients;  they  agree  in  regional  trend  and  magnitude.  The  low- 
order  polynomial  fit  to  the  observations  also  shows  the  same 
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general  trend.  The  older  coefficients  are  selected  for  use  in 
our  present  analysis  in  order  to  be  consistent  with  work  per¬ 
formed  by  the  1AG  Study  Group. 

The  short-wavelength  terrain-predicted  deflections 
(Fig.  C.5-2)  and  the  long-wavelength  deflections  from  the  geo¬ 
potential  model  (Fig.  C.5-3),  when  summed,  provide  an  approxi¬ 
mation  to  the  observed  deflections  of  the  vertical  (Fig.  C.5-1) 
Figure  C.5-4  shows  predicted  deflections  of  the  vertical  which 
are  obtained  by  summing  the  contributions  of  the  terrain  ef¬ 
fects  and  the  geopotential  model.  The  summed  predictions  cor¬ 
relate  well  with  the  direction  of  the  observed  deflections. 
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Predicted  (Geopotential  Model  +  Terrain- 
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F igure  C.5-4 


the  analytic  sciences  corporation 


but  the  magnitudes  are  generally  smaller.  Figure  C. 5-4  also 
agrees  with  "terrain  reduced"  deflections  of  the  vertical  pro¬ 
duced  by  Tscherning  and  Forsberg  (Ref.  20,  Fig.  2,  p.  89). 

Figure  C.5-5  shows  residual  deflections  of  the  ver¬ 
tical  derived  by  subtracting  the  predicted  deflections  (shown 
in  Fig.  C.5-4)  from  the  observed  deflections  of  the  vertical 
(Fig.  C.5-1).  The  residuals  show  systematic  deflection  dif¬ 
ferences  from  4  to  6  arc  sec  in  two  north-south  bands  between 
106  deg  and  106.5  deg  West.  A  possible  source  of  these  sys¬ 
tematic  deflections  is  an  unmodeled  variation  in  density,  which 
could  include,  in  part,  effects  of  isostatic  compensation. 
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Figure  C.5-5  Residual  (Observed-Geopotential  Model  Plus 

Terrain-Predicted)  Deflections  of  the  Vertical 


C-28 


THE  ANALYTIC  SCIENCES  CORPORATION 


Fundamental  statistics  of  the  observed  and  predicted 
deflections  are  tabulated  in  Table  C.5-1.  They  consist  of  the 
rms ,  mean,  and  standard  deviation.  Although  they  are  not  inde 
pendent,  the  rms 


rms 


(C.5-1) 


and  standard  deviation  (o) 


are  useful  in  different  ways,  so  both  are  included. 


TABLE  C.5-1 

STATISTICS  OF  OBSERVED  AND  PREDICTED  DEFLECTIONS 

OF  THE  VERTICAL 


DEFLECTIONS 

1 

FIG. 

NO. 

RMS 

-  _ i 

MEAN 

STANDARD 

DEVIATION 

z 

H 

Z 

u 

Z 

n 

Observed 

5.5-1 

- i 

3.60 

i 

i 

7.78 

l 

-2.44 

1 

-4.56 

2.65 

6.31 

Terrain-Predicted 

5.5-2 

2.20 

4.14 

00 

r-» 

o 

f 

0.12 

2.06 

4.14 

Rapp  '78  Coef. 

5.5-3 

1.26 

4.40 

-1.09 

-4.33 

0.62 

0.77 

Predicted 

fOeopotential  Model 
&  Terrain) 

5.5-4 

2.76 

5.81 

-1.87 

-4.21 

2.03 

4.00 

Residual 
(Observed  -  Rapp 
&  Terrain) 

5.5-5 

j  1.70 

: 

3.15 

-0.57 

-0.35 

1.60 

, 

3.13 
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The  mean  represents  the  very  long  wavelength  portion 
of  each  set  of  deflections,  while  the  standard  deviation  is 
associated  with  the  short  wavelength  amplitude.  The  rms  sup¬ 
plies  a  measure  of  typical  magnitudes  for  the  primary  sets  of 
deflections  and  an  indication  of  goodness  of  fit  in  the  case 
of  deflection  residuals. 


The  success  in  predicting  the  observed  data  is  re¬ 
flected  in  the  statistical  results.  For  instance,  the  mean  of 
the  deflection  residuals  is  nearly  zero  because  the  predic¬ 
tions  account  for  most  of  the  observed  deflections  of  the  ver¬ 
tical.  The  small  rms  of  the  residual  deflections  indicates  a 
high  goodness  of  fit  between  observed  and  modeled  deflections. 
Additional  support  for  the  success  of  the  model  is  shown  by 
the  similarities  between  the  statistical  values  of  the  observed 
and  predicted  deflections. 


The  statistical  analysis  also  provides  a  rough  indi¬ 
cation  of  the  frequency  of  the  terrain-height  and  geopotential 
model  components  of  the  predicted  deflections.  The  near-zero 
mean  of  the  terrain-predicted  deflections  show  that  their  con¬ 
tribution  is  primarily  short  wavelength.  The  low  standard 
deviation  indicates  a  long  wavelength  contribution. 


The  coefficient  of  correlation,  (r),  is  computed  via: 


(C.5-3) 


to  quantify  the  quality  of  the  linear  fit  between  the  observed 
and  predicted  deflections.  For  the  east  component,  r  =  0.91; 
for  the  north  component,  r  =  0.80.  These  high  positive  values 
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clearly  show  how  well  the  predicted  deflections  approximate 
the  observed  deflection  values.  Interestingly,  most  of  the 
high  correlation  is  due  to  the  terrain-predicted  values  be¬ 
cause  the  coefficient  of  correlation  between  the  observed  de¬ 
flections  and  deflections  predicted  from  terrain  only  is  vir¬ 
tually  identical  to  the  coefficient  values  listed  above. 


A  more  detailed  view  of  the  correlation  between  ob¬ 
served  and  modeled  deflections  is  provided  by  the  scatter  dia¬ 
grams  shown  in  Figs.  C.5-6  and  C.5-7.  The  fact  that  the  least- 
squares  regression  line  for  the  east  vertical  deflections  has 
a  slope  greater  than  one  displays  the  systematic  under-estimation 
of  the  east  component  of  the  observed  deflections  by  the  modeled 
deflections,  as  shown  in  Fig.  C.5-5.  If  isostatic  compensa¬ 
tion  were  added,  the  underestimation  might  become  even  more 
pronounced.  Therefore,  the  result  points  to  systematic  lateral 
density  variations. 
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Figure  C.5-6  Scatter  Diagram  of  North  Vertical 
Deflection  Component 


C-31 


THE  ANALYTIC  SCIENCES  CORPORATION 


A-4Q243 


-36.0  -30.0  -26.0  -20.0  -16.0  -10.0  -5.0  0.0  5.0  10.0  16.0 

PREDICTED  DEFLECTIONS  OF  THE  VERTICAL  (arc  sacs) 


Figure  C.5-7  Scatter  Diagram  of  East  Vertical 
Deflection  Component 


In  conclusion,  it  may  be  stated  that  observed  deflec¬ 
tions  of  the  vertical  can  be  approximated  reasonably  well  using 
deflections  predicted  from  local  terrain  and  a  global  geopoten¬ 
tial  model.  Of  course,  with  gravimetric  deflection  predictions, 
far  more  accuracy  is  attainable,  as  shown  in  the  IAG  Study 
Group  reports  (Refs.  20  and  21).  The  work  reported  here  is 
complementary  to  the  study  group  research  in  providing  a  more 
thorough  examination  of  the  terrain  contribution  and  the  novel 
examination  of  correlation  and  scatter  diagrams. 
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